{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4ed546c8-cbad-4721-ba6b-7d185222990e",
   "metadata": {},
   "source": [
    "# Example 4\n",
    "No source, Multiply-connected domain, Dirichlet/Neumann problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "228a690d-545c-48b5-a2d5-6a3e40d7963b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import calfem.geometry as cfg\n",
    "import calfem.mesh as cfm\n",
    "import calfem.vis as cfv\n",
    "import numpy as np\n",
    "from toolkits import gauss_line_int, gauss_tri_int\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b47a234-100f-4115-a5ee-670e2b590fdd",
   "metadata": {},
   "source": [
    "## 1. Geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "1002c549-ae11-4ae3-b712-f909f8ba093c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of holes\n",
    "nH = 1\n",
    "# Mesh size factor\n",
    "factor = 20\n",
    "\n",
    "# externel boundary ID\n",
    "dir_id_1 = 100\n",
    "dir_id_2 = 200\n",
    "dir_id_3 = 300\n",
    "neu_id_1 = 400\n",
    "# hole boundaryd ID\n",
    "hole_neu_id_1 = 500\n",
    "hole_neu_id_2 = 600\n",
    "hole_neu_id_3 = 700\n",
    "hole_neu_id_4 = 800\n",
    "\n",
    "# Geometry definition\n",
    "g = cfg.Geometry()\n",
    "g.point([0.0, 0.0], ID=0)\n",
    "g.point([1.0, 0.0], ID=1)\n",
    "g.point([1.0, 1.0], ID=2)\n",
    "g.point([0.0, 1.0], ID=3)\n",
    "g.point([0.4, 0.4], ID=4)\n",
    "g.point([0.6, 0.4], ID=5)\n",
    "g.point([0.6, 0.6], ID=6)\n",
    "g.point([0.4, 0.6], ID=7)\n",
    "g.spline([0, 1], el_on_curve=factor, marker=dir_id_1, ID=0)\n",
    "g.spline([1, 2], el_on_curve=factor, marker=neu_id_1, ID=1) \n",
    "g.spline([2, 3], el_on_curve=factor, marker=dir_id_2, ID=2) \n",
    "g.spline([3, 0], el_on_curve=factor, marker=dir_id_3, ID=3)\n",
    "g.spline([4, 5], el_on_curve=factor/4, marker=hole_neu_id_1, ID=4)\n",
    "g.spline([5, 6], el_on_curve=factor/4, marker=hole_neu_id_2, ID=5)\n",
    "g.spline([6, 7], el_on_curve=factor/4, marker=hole_neu_id_3, ID=6)\n",
    "g.spline([7, 4], el_on_curve=factor/4, marker=hole_neu_id_4, ID=7)\n",
    "g.surface([0, 1, 2, 3],[[4,5,6,7]])\n",
    "cfv.drawGeometry(g)\n",
    "cfv.showAndWait()\n",
    "\n",
    "# g.point([0.0, 0.0], ID=0)\n",
    "# g.point([1.0, 0.0], ID=1)\n",
    "# g.point([1.0, 1.0], ID=2)\n",
    "# g.point([0.0, 1.0], ID=3)\n",
    "# g.point([0.5, 0.5], ID=4)\n",
    "# g.point([0.6, 0.5], ID=5)\n",
    "# g.point([0.5, 0.6], ID=6)\n",
    "# g.point([0.4, 0.5], ID=7)\n",
    "# g.point([0.5, 0.4], ID=8)\n",
    "# g.spline([0, 1], el_on_curve=factor, marker=dir_id_1, ID=0)\n",
    "# g.spline([1, 2], el_on_curve=factor, marker=neu_id_1, ID=1) \n",
    "# g.spline([2, 3], el_on_curve=factor, marker=dir_id_2, ID=2) \n",
    "# g.spline([3, 0], el_on_curve=factor, marker=dir_id_3, ID=3)\n",
    "# g.circle([5, 4, 6], el_on_curve=factor/4, ID=4)\n",
    "# g.circle([6, 4, 7], el_on_curve=factor/4, ID=5)\n",
    "# g.circle([7, 4, 8], el_on_curve=factor/4, ID=6)\n",
    "# g.circle([8, 4, 5], el_on_curve=factor/4, ID=7)\n",
    "# g.surface([0, 1, 2, 3],[[4,5,6,7]])\n",
    "# # g.structuredSurface([0, 1, 2, 3])\n",
    "# Display mesh\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90c4263d-6fe7-484d-8ca4-8a590e25a3b6",
   "metadata": {},
   "source": [
    "## 2. Meshing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "07931634-1ddc-4c54-9b23-e70b883c5d48",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : GMSH -> Python-module\n",
      "[[ 4 99]\n",
      " [99 98]\n",
      " [98 97]\n",
      " [97 96]\n",
      " [96  7]\n",
      " [ 7 95]\n",
      " [95 94]\n",
      " [94 93]\n",
      " [93 92]\n",
      " [92  6]\n",
      " [ 6 91]\n",
      " [91 90]\n",
      " [90 89]\n",
      " [89 88]\n",
      " [88  5]\n",
      " [ 5 87]\n",
      " [87 86]\n",
      " [86 85]\n",
      " [85 84]\n",
      " [84  4]]\n",
      "[[ 5 87]\n",
      " [87 86]\n",
      " [86 85]\n",
      " [85 84]\n",
      " [84  4]]\n",
      "[[ 6 91]\n",
      " [91 90]\n",
      " [90 89]\n",
      " [89 88]\n",
      " [88  5]]\n",
      "[[ 7 95]\n",
      " [95 94]\n",
      " [94 93]\n",
      " [93 92]\n",
      " [92  6]]\n",
      "[[ 4 99]\n",
      " [99 98]\n",
      " [98 97]\n",
      " [97 96]\n",
      " [96  7]]\n",
      "Dirichlet !\n"
     ]
    }
   ],
   "source": [
    "mesh = cfm.GmshMesh(g,return_boundary_elements=True)\n",
    "mesh.elType = 2          # Degrees of freedom per node.\n",
    "mesh.dofsPerNode = 1     # Factor that changes element sizes.\n",
    "# mesh.elSizeFactor = 0.5 # Element size Factor\n",
    "\n",
    "coords, edof, dofs, bdofs, elementmarkers, boundaryElements  = mesh.create()\n",
    "\n",
    "#############################\n",
    "#                           #\n",
    "# Elements, Edges and Nodes #\n",
    "#                           #\n",
    "#############################\n",
    "element = edof - 1\n",
    "node = coords\n",
    "edge = dict()\n",
    "for i, e in enumerate(element):\n",
    "    edge[(e[0], e[1])] = i\n",
    "    edge[(e[1], e[2])] = i\n",
    "    edge[(e[2], e[0])] = i\n",
    "\n",
    "######################\n",
    "#                    #\n",
    "#   Externel edges   #\n",
    "#                    #\n",
    "######################\n",
    "dir_edge_1 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_1]])-1\n",
    "dir_edge_2 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_2]])-1\n",
    "dir_edge_3 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_3]])-1\n",
    "\n",
    "neu_edge_1 = np.array([edge['node-number-list'] for edge in boundaryElements[neu_id_1]])-1\n",
    "\n",
    "dir_ids = [dir_id_1, dir_id_2, dir_id_3]\n",
    "neu_ids = [neu_id_1]\n",
    "\n",
    "dir_edges = set([tuple(np.array(edge['node-number-list'])-1) for marker in dir_ids for edge in boundaryElements[marker]])\n",
    "neu_edges = set([tuple(np.array(edge['node-number-list'])-1) for marker in neu_ids for edge in boundaryElements[marker]])\n",
    "\n",
    "######################\n",
    "#                    #\n",
    "#     Hole edges     #\n",
    "#                    #\n",
    "######################\n",
    "# first hole\n",
    "hole_neu_edge_1 = np.array([edge['node-number-list'] for edge in boundaryElements[hole_neu_id_1]])-1\n",
    "hole_neu_edge_2 = np.array([edge['node-number-list'] for edge in boundaryElements[hole_neu_id_2]])-1\n",
    "hole_neu_edge_3 = np.array([edge['node-number-list'] for edge in boundaryElements[hole_neu_id_3]])-1\n",
    "hole_neu_edge_4 = np.array([edge['node-number-list'] for edge in boundaryElements[hole_neu_id_4]])-1\n",
    "\n",
    "hole_neu_edges = [hole_neu_edge_1, hole_neu_edge_2, hole_neu_edge_3, hole_neu_edge_4]\n",
    "\n",
    "ccw_hole_neu_edges = np.array([edge[::-1] for edges in hole_neu_edges[::-1] for edge in edges[::-1]])\n",
    "\n",
    "ccw_hole_neu_edge_1 = np.array([edge[::-1] for edge in hole_neu_edge_1[::-1]])\n",
    "ccw_hole_neu_edge_2 = np.array([edge[::-1] for edge in hole_neu_edge_2[::-1]])\n",
    "ccw_hole_neu_edge_3 = np.array([edge[::-1] for edge in hole_neu_edge_3[::-1]])\n",
    "ccw_hole_neu_edge_4 = np.array([edge[::-1] for edge in hole_neu_edge_4[::-1]])\n",
    "\n",
    "print(ccw_hole_neu_edges)\n",
    "print(ccw_hole_neu_edge_1)\n",
    "print(ccw_hole_neu_edge_2)\n",
    "print(ccw_hole_neu_edge_3)\n",
    "print(ccw_hole_neu_edge_4)\n",
    "\n",
    "#################################\n",
    "#                               #\n",
    "#  Crossing elements and edges  #\n",
    "#                               #\n",
    "#################################\n",
    "direction = np.array([0, -1])\n",
    "Bh = []\n",
    "Eh = {}\n",
    "\n",
    "Bh.append(tuple(hole_neu_edge_1[0]))\n",
    "\n",
    "# crossing flag \n",
    "connected = False\n",
    "\n",
    "######################################################\n",
    "# TODO: New structures of Eh, Bh for multiple holes. #\n",
    "######################################################\n",
    "while not connected:\n",
    "    cross_edge = Bh[-1]\n",
    "    # print(cross_edge)\n",
    "    # print(node[cross_edge[0]])\n",
    "    # print(node[cross_edge[1]])\n",
    "    if (cross_edge not in dir_edges) and (cross_edge not in neu_edges):\n",
    "        cross_edge = cross_edge[::-1]\n",
    "        ei = edge[cross_edge]\n",
    "        \n",
    "        i, j, k = element[ei]\n",
    "        order = [i,j,k]  \n",
    "        \n",
    "        local_edges = set([(i,j), (j,k), (k,i)])\n",
    "        local_edges.remove(cross_edge)\n",
    "        e1 = np.array(local_edges.pop())\n",
    "        e2 = np.array(local_edges.pop())\n",
    "        \n",
    "        e1n = node[e1[1]] - node[e1[0]]\n",
    "        e1n = e1n/np.linalg.norm(e1n)\n",
    "        \n",
    "        e2n = node[e2[1]] - node[e2[0]]\n",
    "        e2n = e2n/np.linalg.norm(e2n)\n",
    "        \n",
    "        d1 = np.abs(np.dot(direction, e1n))\n",
    "        d2 = np.abs(np.dot(direction, e2n))\n",
    "        if d1 < d2:\n",
    "            ae = np.zeros(3)\n",
    "            L1 = np.linalg.norm(node[cross_edge[0]] - node[cross_edge[1]])\n",
    "            L2 = np.linalg.norm(node[e1[0]] - node[e1[1]])\n",
    "            ae[order.index(cross_edge[0])]= 1/L1\n",
    "            ae[order.index(e1[0])]= -1/L2\n",
    "            Eh[ei] = ae\n",
    "            Bh.append(tuple(e1))\n",
    "        else:\n",
    "            ae = np.zeros(3)\n",
    "            L1 = np.linalg.norm(node[cross_edge[0]] - node[cross_edge[1]])\n",
    "            L2 = np.linalg.norm(node[e2[0]] - node[e2[1]])\n",
    "            ae[order.index(cross_edge[0])]= 1/L1\n",
    "            ae[order.index(e2[0])]= -1/L2\n",
    "            Eh[ei] = ae\n",
    "            Bh.append(tuple(e2))\n",
    "        # print('Element={}, ae={}, next <{},{}>'.format(ei, ae, *Bh[-1]))            \n",
    "        # print(e1,'\\t',e1n)\n",
    "        # print(e2,'\\t',e2n)\n",
    "        # print(ae)\n",
    "        # print(d1,d2)\n",
    "        # print(\"----------------------------------\")\n",
    "        \n",
    "            \n",
    "    else:\n",
    "        if cross_edge in neu_edges:\n",
    "            print('Neumann ! ')\n",
    "        else:\n",
    "            print('Dirichlet !')\n",
    "        connected = True\n",
    "    \n",
    "# print(Eh)\n",
    "        \n",
    "# np.all()\n",
    "\n",
    "######################\n",
    "#                    #\n",
    "#  Numbers of sets   #\n",
    "#                    #\n",
    "######################\n",
    "E = len(element)\n",
    "N = len(node)\n",
    "num_dof = N + nH\n",
    "\n",
    "# Display mesh\n",
    "cfv.draw_mesh(coords=coords, edof=edof, dofs_per_node=mesh.dofs_per_node,\n",
    "              el_type=mesh.el_type,filled=True,title=\"Example 2 - Mesh\")\n",
    "cfv.showAndWait()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edf91826-efaf-416e-9306-d783aa7ea2c6",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. Boundary condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "cb26e870-b29d-41b5-bb21-0c83e022121b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Externel boundary condition\n",
    "# ########## Dirichlet ##########\n",
    "# def dir1(x, y):\n",
    "#     return np.zeros_like(x)\n",
    "\n",
    "# def dir2(x, y):\n",
    "#     return np.sin(x)/np.sin(1)\n",
    "\n",
    "# def dir3(x, y):\n",
    "#     return np.zeros_like(x)\n",
    "# ########## Neumann ############\n",
    "# def neu1(x, y):\n",
    "#     return (1/np.tan(1))*np.sinh(y)/np.sinh(1)\n",
    "\n",
    "# # Hole's boundary condition\n",
    "# ########## Dirichlet ##########\n",
    "# # ...\n",
    "# ########## Neumann ############\n",
    "# def hole_neu1(x, y):\n",
    "#     return (np.sin(x)/np.sin(1))*(np.cosh(0.4)/np.sinh(1))\n",
    "\n",
    "# def hole_neu2(x, y):\n",
    "#     return -(np.cos(0.6)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "# def hole_neu3(x, y):\n",
    "#     return -(np.sin(x)/np.sin(1))*(np.cosh(0.6)/np.sinh(1))\n",
    "\n",
    "# def hole_neu4(x, y):\n",
    "#     return (np.cos(0.4)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "\n",
    "\n",
    "# Externel boundary condition\n",
    "########## Dirichlet ##########\n",
    "def dir1(x, y):\n",
    "    return exact_u(x,y)\n",
    "\n",
    "def dir2(x, y):\n",
    "    return exact_u(x,y)\n",
    "\n",
    "def dir3(x, y):\n",
    "    return exact_u(x,y)\n",
    "########## Neumann ############\n",
    "def neu1(x, y):\n",
    "    return exact_ux(x,y)\n",
    "\n",
    "# Hole's boundary condition\n",
    "########## Dirichlet ##########\n",
    "# ...\n",
    "########## Neumann ############\n",
    "def hole_neu1(x, y):\n",
    "    return exact_uy(x,y)\n",
    "\n",
    "def hole_neu2(x, y):\n",
    "    return -exact_ux(x,y)\n",
    "\n",
    "def hole_neu3(x, y):\n",
    "    return -exact_uy(x,y)\n",
    "\n",
    "def hole_neu4(x, y):\n",
    "    return exact_ux(x,y)\n",
    "\n",
    "def exact_u(x,y):\n",
    "    return x*(1-x)+y*(y-1)\n",
    "\n",
    "def exact_ux(x,y):\n",
    "    return 1-2*x\n",
    "\n",
    "def exact_uy(x,y):\n",
    "    return 2*y-1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5fe1912-f8d3-4895-9e96-3156355eb71c",
   "metadata": {},
   "source": [
    "## 3. Assemble flexibility matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "fb18dbeb-89fc-4566-b9bc-7627f4bef33f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "21 [155 379 516  -1]\n",
      "40 [123 367 368  -1]\n",
      "102 [158 382 463  -1]\n",
      "109 [382 367 463  -1]\n",
      "142 [368 367 382  -1]\n",
      "153 [379 329 516  -1]\n",
      "172 [516 329 584  -1]\n",
      "380 [149 154 155  -1]\n",
      "448 [155 154 156  -1]\n",
      "487 [ 84   4 113  -1]\n",
      "582 [156 154 158  -1]\n",
      "612 [123 368 552  -1]\n",
      "664 [113   4 140  -1]\n",
      "715 [113 140 329  -1]\n",
      "728 [156 158 463  -1]\n",
      "810 [149 155 516  -1]\n",
      "839 [ 15 123 552  -1]\n",
      "929 [329 140 584  -1]\n",
      "1016 [ 14  15 552  -1]\n",
      "size: 597 x 597\n",
      "rank: 596\n"
     ]
    }
   ],
   "source": [
    "F = np.zeros((num_dof, num_dof))\n",
    "\n",
    "for i, e in enumerate(element):\n",
    "    # print(e)\n",
    "    # print(node.shape)\n",
    "    x1,y1 = node[e][0]\n",
    "    x2,y2 = node[e][1]\n",
    "    x3,y3 = node[e][2]\n",
    "    \n",
    "    A_e = 1/2*np.linalg.det(np.c_[node[e], np.ones((3,1))])\n",
    "    \n",
    "    Phi_e = -(1/(2*A_e))*np.array([[x2-x3,x3-x1,x1-x2],\n",
    "                                  [y2-y3,y3-y1,y1-y2]])\n",
    "    \n",
    "    \n",
    "    if i in Eh:\n",
    "        L1 = np.sqrt((x2-x1)**2+(y2-y1)**2)\n",
    "        L2 = np.sqrt((x3-x2)**2+(y3-y2)**2)\n",
    "\n",
    "        left = np.linalg.inv(-np.array([[(y1-y2)/L1,(x2-x1)/L1],\n",
    "                              [(y2-y3)/L2,(x3-x2)/L2]]))\n",
    "\n",
    "        Phi_L = np.c_[left, np.zeros((2,1))]\n",
    "        \n",
    "        # print(Eh[i])\n",
    "        Phi_e = np.c_[Phi_e, Phi_L @ Eh[i]]\n",
    "        F_e = A_e * Phi_e.T @ Phi_e\n",
    "        \n",
    "        ### Todo: multiple holes ###\n",
    "        ey = np.append(e, -1)\n",
    "\n",
    "        F[np.ix_(ey,ey)] = F[np.ix_(ey,ey)] + F_e\n",
    "        print(i, ey)\n",
    "    else:\n",
    "        F_e = A_e * Phi_e.T @ Phi_e\n",
    "        \n",
    "        F[np.ix_(e,e)] = F[np.ix_(e,e)] + F_e\n",
    "\n",
    "FFF = F.copy()\n",
    "print('size: {} x {}'.format(*F.shape))\n",
    "print('rank: {}'.format(np.linalg.matrix_rank(F)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a01a3f7-c7e0-40b2-aeef-6ac3d2787a52",
   "metadata": {},
   "source": [
    "## 4. Assemble kinematic vector (with source, Dirichlet problem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "c855dbeb-fe00-4b49-8f23-8514e2fd3e53",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 2]\n",
      "[99, 98, 97, 96, 7, 95, 94, 93, 92, 6, 91, 90, 89, 88, 5, 87, 86, 85, 84]\n",
      "(4, 99)\n",
      "(99, 98)\n",
      "(98, 97)\n",
      "(97, 96)\n",
      "(96, 7)\n",
      "(7, 95)\n",
      "(95, 94)\n",
      "(94, 93)\n",
      "(93, 92)\n",
      "(92, 6)\n",
      "(6, 91)\n",
      "(91, 90)\n",
      "(90, 89)\n",
      "(89, 88)\n",
      "(88, 5)\n",
      "(5, 87)\n",
      "(87, 86)\n",
      "(86, 85)\n",
      "(85, 84)\n"
     ]
    }
   ],
   "source": [
    "du = np.zeros(num_dof)\n",
    "\n",
    "##########################\n",
    "#   Dirichlet boundary   #\n",
    "##########################\n",
    "for edge in dir_edge_1:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir1, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "\n",
    "for edge in dir_edge_2:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir2, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "\n",
    "for edge in dir_edge_3:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir3, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "\n",
    "##########################################\n",
    "#    Neumann boundary (without source)   #\n",
    "##########################################\n",
    "s_N = np.zeros(num_dof)\n",
    "A = np.eye(num_dof)\n",
    "\n",
    "num_red = len(neu_edge_1) + len(ccw_hole_neu_edges)\n",
    "num_eff = num_dof - num_red\n",
    "\n",
    "# 1. Externel edge\n",
    "# effective variable (s1) is at first node of first Neumann edge \n",
    "# (edge's nodes are ordered counterclockwise within the local element).\n",
    "ext_eff_dofs = [neu_edge_1[0][0]]\n",
    "ext_red_dofs = [edge[1] for edge in neu_edge_1] \n",
    "\n",
    "# 2. Hole edge\n",
    "# last edge should be excluded for closed Neumann boundary. \n",
    "hole_eff_dofs = [hole_neu_edge_1[0][0]]\n",
    "hole_red_dofs = [edge[1] for edge in ccw_hole_neu_edges[:-1]]\n",
    "\n",
    "red_dofs = ext_red_dofs + hole_red_dofs\n",
    "print(ext_red_dofs)\n",
    "print(hole_red_dofs)\n",
    "\n",
    "\n",
    "A[np.ix_(ext_red_dofs, ext_eff_dofs)] = 1\n",
    "A[np.ix_(hole_red_dofs, hole_eff_dofs)] = 1\n",
    "\n",
    "A = np.delete(A, red_dofs, axis=1)\n",
    "\n",
    "# Neumann terms\n",
    "####### TODO : crossing edge #########\n",
    "# externel (counter-clockwise)\n",
    "fl_bar = 0\n",
    "for edge in neu_edge_1:\n",
    "    i, j = edge\n",
    "    fl_bar = fl_bar + gauss_line_int(neu1, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "\n",
    "# hole (clockwise)\n",
    "# last edge should be excluded for closed Neumann boundary. \n",
    "fl_bar = 0\n",
    "for edge in ccw_hole_neu_edge_4:\n",
    "    i, j = edge\n",
    "    print((i,j))\n",
    "    fl_bar = fl_bar + gauss_line_int(hole_neu4, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "    \n",
    "for edge in ccw_hole_neu_edge_3:\n",
    "    i, j = edge\n",
    "    print((i,j))\n",
    "    \n",
    "    fl_bar = fl_bar + gauss_line_int(hole_neu3, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "    \n",
    "for edge in ccw_hole_neu_edge_2:\n",
    "    i, j = edge\n",
    "    print((i,j))\n",
    "    \n",
    "    fl_bar = fl_bar + gauss_line_int(hole_neu2, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "    \n",
    "for edge in ccw_hole_neu_edge_1[:-1]:\n",
    "    i, j = edge\n",
    "    print((i,j))\n",
    "    \n",
    "    fl_bar = fl_bar + gauss_line_int(hole_neu1, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "    \n",
    "# Source \n",
    "# Particular solution q_0\n",
    "# def qx_0(x, y):\n",
    "#     return x\n",
    "\n",
    "# def qy_0(x, y):\n",
    "#     return np.zeros_like(x)\n",
    "\n",
    "# for i, e in enumerate(element):\n",
    "#     x1,y1 = node[e][0]\n",
    "#     x2,y2 = node[e][1]\n",
    "#     x3,y3 = node[e][2]\n",
    "    \n",
    "#     A_e = 1/2*np.linalg.det(np.c_[node[e], np.ones((3,1))])\n",
    "    \n",
    "#     Phi_e = (1/(2*A_e))*np.array([[x2-x3,x3-x1,x1-x2],\n",
    "#                                   [y2-y3,y3-y1,y1-y2]])\n",
    "    \n",
    "#     int_qx = gauss_tri_int(qx_0, x1, y1, x2, y2, x3, y3, degree=2)\n",
    "#     int_qy = gauss_tri_int(qy_0, x1, y1, x2, y2, x3, y3, degree=2)\n",
    "    \n",
    "#     par = Phi_e.T @ np.array([int_qx, int_qy])\n",
    "    \n",
    "#     du[e] = du[e] - par\n",
    "\n",
    "# Neumann boundary with source.\n",
    "# ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5619f835-d960-4905-a43c-5b77ff0dcc40",
   "metadata": {},
   "source": [
    "## 5. Constraint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "0dba5f60-fa70-4d41-ac50-72d386d3aa69",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.9081958235744878e-17\n"
     ]
    }
   ],
   "source": [
    "FU = F.copy()\n",
    "dU = du.copy()\n",
    "\n",
    "F = A.T @ FU @ A\n",
    "d = A.T @ (dU - FU@s_N)\n",
    "\n",
    "r1 = 0\n",
    "\n",
    "for edge in hole_neu_edge_1:\n",
    "    i, j = edge\n",
    "    r1 = r1 + gauss_line_int(hole_neu1, *node[i], *node[j])\n",
    "\n",
    "for edge in hole_neu_edge_2:\n",
    "    i, j = edge\n",
    "    r1 = r1 + gauss_line_int(hole_neu2, *node[i], *node[j])\n",
    "\n",
    "for edge in hole_neu_edge_3:\n",
    "    i, j = edge\n",
    "    r1 = r1 + gauss_line_int(hole_neu3, *node[i], *node[j])\n",
    "\n",
    "for edge in hole_neu_edge_4:\n",
    "    i, j = edge\n",
    "    r1 = r1 + gauss_line_int(hole_neu4, *node[i], *node[j])\n",
    "    \n",
    "print(r1)\n",
    "\n",
    "F[-1,:] = 0\n",
    "F[-1,-1] = 1\n",
    "d[-1] = r1\n",
    "\n",
    "F[0,:] = 0\n",
    "F[0,0] = 1\n",
    "d[0] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2c04f820-d63c-4dbd-9f7d-00b07d8bec82",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x17ce59dcac0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkcAAAJCCAYAAADKjmNEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAB360lEQVR4nO3dX4wlWX4X+O+JiHvzb1VmdVVlTlZmdXW3PRKah8VYLe8g/ABGsLYXMX6wLLOsPLJGmhcejGDFGl5YVot2ecFgsbJkYYsxAoxlYD1CaMVobAlebDyDjW1sWDfdXVVZVZ1ZVZ1/qvLPvfHn7EPeczvy1v0Tf86JOOfE9yO1OvNmVt64ESfO+Z3f+RNCSgkiIiIiuhK0fQBERERENmFwRERERJTD4IiIiIgoh8ERERERUQ6DIyIiIqIcBkdEREREOUaCIyHE9wsh/qsQ4gMhxE+aeA8iIiIiE4TufY6EECGA/w/AnwGwD+A3AfwFKeXva30jIiIiIgNMZI6+B8AHUsoPpZRDAL8I4EsG3oeIiIhIu8jA39wF8Dj3/T6A/37ePwiCQAYBpz8RERFRc9I0fSGlvDv5uongqBAhxFcBfHX0NR49eoQvfOEL2v6+lBJSStgadMVxjCAIEIYhAODy8hL9fr/w8Z6fn2N1dXX8/eXlJZaXl6/9jpQSg8EAy8vLyLIMcRxjaWkJaZoiTVP0+319H4i8MK0ckb8m65Ei8vVHkiSQUqLX6137neFwiDAMx/Wb6WPyyXA4RBRFtdsuHedxMBggCII3ru+kV69eYX19HUKIa6+fnZ1hbW1t6vdSSrx+/Ro3btyodYx1nZycPJz2uok5R38cwP8mpfwfRt//dQCQUv6fs/5NFEVyfX193HAnSaL1mOrKsgwArA20iIiIbJamaaVg2bSTk5NvSynfn3zdROboNwF8XgjxLoAnAH4UwP9U5B+maYoPP/wQDx48gE0PxG0yKLq4uBgHY8BVVk1F41JKrT2qNE0xHA7H36segotBYBzHkFI6mQ0bDodOHjeRKRcXF+M2QAiBlZWVlo/oTRcXFwBg5bHZqExglCQJzs7O0Ov1jGYR52XEtGeOAEAI8YMA/h6AEMDPSyn/9rzfD8NQqtSaEAI3b97E0dGRk400EREtlmUZ63iPuTJEPytzZCQ4KksNqylZlmFzcxOnp6ctHhUREZmSJAmiqLVpr50lpUSSJAvnEXXFrODIyrA9CAKcnp7ixo0bVg2vuUpKiTRN2z4MIqIxBkZkMyuDI0UFSFRPft4SERG1S634a4MQglmjAqwOjoQQeP369bWlgL5RS2FN49g+EZEdgiBgnWw5J67O2dmZt/teMKNDRK7LsuzaKluaz+dsvi9TYZwIjoCrDa3W1ta8i7bDMNR2k/hSKInILT439lTOYDBo+xAKm9dmOhVpXFxcYHl5mTfhDJx0TURtYHD0ma53Ul1KYOT3+ZvkzqfAVepWBUguXQBTJm9CVk5ERO3qeifVpQ1tl5aWZv7MuQhDPS+szHPIfDX5mBUbt2YnInd0PeuhAzupfnAyusiybPzIhS4HSHU/OytCIsqz7bmWLmIn1Q/ORhYqQHL1WWA61LkJubqEiCYFQcBOExEcDo6AqwY+jmNnA6Q2KyFOoCSiSWEYMntEBMeDI8DtACnLstYCJCGEc+eLiMzL1wvMLlNXedE6qgApiiKnxnubCk5YwRFRUfk6lHUHdZUXwRFwdROnaYooipzJiDQxtCWlZAVHRJW4UpdOklJy7hTV4mbJnyFNUyczSCZx+IzKYjBNiqt1B4Mjqsuakq+rIGdZhiRJGCDluFrBUTsYHJHruODETi5N9rem1dQZ5asAqev7IBFVwU4F6dLmghMGR+2L4/ja9y5dE2siB91BTJqmGAwGTq5iI2qTSxUY2a3NFbk0n5TS+KNOJttelzpeXkcNahVbVzNIHB4hojYFQcBg22Km2wiXgqFJ3kcM+UeNdO0mZY+NiNrUtTq3jLafUiCEQBRFrb2/7bwPjoCrQjgYDLC8vGz1zZokidaAxuWonYiIzLK5PWxbJ4Ij4CqLcnl5iZWVlbYPhYjIO5xfVE4QBJ2c7uGKTl0ZKSXOz8+xtrZmZaGMooiRPBEROW0wGLR9CLXZFyE04OLiAqurq20fhjamVxwQES3Cydek+JBB7GRwlGUZXr9+jbW1tbYPRQsfCiIRkW+aWC5vo6WlpbYPobZOBkfK2dkZ1tfXnQsuJncZ5YoDIiI7uda+6OBDBrHTwREAvHr1Cjdv3nSqALt0rETkhsndjKk+Lpd3V+eDIyEETk9PnQqQer1e5X/b1TQvERFRUZ0PjoDPAqSNjQ1nAiQiorLmZYfqdLqIfMPgaEQIgZOTE2xsbBjdtdSGXVG5OSQREdFsDI5yhBA4Pj7GrVu3+FwyIvJOPjvEOUZEszE4mhAEAY6OjowFSE3sispKj4gW8WFFEZEpDI6myAdInLxMRD7ybRUVF5uQTgyOZgiCAJ9++ineeust54bYOLGSyprcO4uIqMusCo5sC0LCMBxnkIiIyF5cbGI31563ZlVwZKMgCBpZxUbUJt+GWIhclGWZt9vJuPa5rAqOTE9UruPo6Aibm5vOXWAiInKDrvbFxmHy5eXltg+hFHujEcvkM0guBEg23hxERDRbGIZaVhFylKM+BkclqI0iXXrUCBFRVVJKdrQc1O/32z4E5zE4Kin/LDZdkiTRHmxxDgkR6cD9kKiLGBxVIITAq1evcOPGDS0VBysfIrIRV4CRKbavXmNwVMPr16+xtrZWO7jRNc5MRESkQ5qmRjfVtHkBFsDgqBYpJc7OzrC6utraMXA+ABHNw/mRVJXJTrvtmxUzOKopHyA1kf2ZrOiYcSKiediBoirCMLQ+u2NSdz+5ZhcXF1heXjYerExWdJwPQETzuNTAMctFtnDnrrGclBKXl5fGA6S6FR0rH6JucakDxSwX2YLBkUYqQFpaWjLWW6tT0UkpuTkYEVnLpSwXXfG1TWFJ1ExKieFwiKWlpZmBTJvZG85RIiJbuZTloism9umzAYMjA7Isw3A4RBRFU3tCbWVwhBDsmRFRZT42glRPEARedrrZUhqSpiniOEav13sjIBFCGCtMvqY4iahdUkqj+96Qm3x9GgODI4OyLJsaIJkMjlh5EZEpOjLP7MCRCxgcGZYPkJoYT+eYPVXFRovm0TUsz3JGLrAiOPJ9HDvLMiRJgiiKjAcvnFNEVbHRoiawjuoG17dl8HOwsAbVQOi+gdVw18rKCuI4tv6he9Q9zDr6T9VvUsrWFmgwOOoG1ydpW1FKXT+J06gVafmsmJqkvbKygqWlpRaPjuhNPt6HtrAtO171Wk/WadQdaZpOvfazMs6ud7asyRypeTlt09WrEUJMLUgqY7SysnLteyLylwpGVMamLXXrt8ljl1KOPxODa78JIZCmKcIwfONaq9d9Yk1wVOemVVka2/ZbmPWZVEC0uro63jSSiPxnU/00i5oCIKVEEARvrLTNU9+78LmoHtW+Tl7rIAi8zCZaMawGzE/BJUkyd4m6i5sbDgYDnJ2d4caNG+j3+20fDhERAIwDomkZgmlUltzHBtKULMucXAAxqzwIIbzr5DsRUYRhuDD4aTqtuyhgK2I4HOLVq1cMkIjIGqouVf+pem5eY85htXJsG+XQoexmkHWDwzRNjQaY1gRHl5eXM39m441XJGArYjgc4vT0FDdv3rRizhURUZ6q59T/84tNmC2qzrY2ra6y7WEQBLWCG/XvTQVI1gRHTWROkiTRtveCzoAtjmOcnJxgc3OTARIRWWXaPKNp808YLJXj6tDaPGUXGNVJMAghtCUpprEmOFIf8OLiwth7hGFo7Yz6OI5xfHyMzc1Nb59VYxIr5XaZemyN6xvJ6eBK0GFjht9mk5PdfaAryVH04ewmy5t1V8ZkJG37zRvHMY6OjnDnzp3WgrhZe1nYzubr2gWmyis7CvbXW10Wx7HX71eWrnJqwyIr64KjrkuSBM+fP8fW1lYrAVLRFSo2klIy00BEjWhjb74oirxbFWYrBkcWStMUh4eH2Nraaj16dokQwslMg4uZuml0D60x0P2MK0NrXdLG/FAhRGdWNhcdWjPFupbX1ayFbmma4uDgANvb2wyQCnC5IXWpzM9L6+vOdLoU6Jouf10bWnP5fiY92h5as67V7VIFsEiWZTg4OMDOzk4j58X0vhEmudSQuszl1ZRSSmMTx5sofz6ubppFnU/b59iQ2+YtALMuOKLrsizD06dPsbu7azxAMrkssgnsbdI8rj//adrqJlcXUBRhy/M2yV/qGafTWNcS+nqj1yGlxJMnT7C3t8fM2hyuZo9c6x27dryK6fLRRnDu8gKKRVRgxE4PtcG64Gh1dbXtQ6hESok4jo2l7aWU2N/fx97enpG/T+1xrXfs2vE2xdbgPEkSZwLaacfZ1HllEEZ51gVHTdNVcQgh0Ov1jKbtpZR4/PgxHjx4YOw9iMgvURQVCmhNdezKaDPwNhGEmZznpkuWZdYEhirJYANrgqPz8/NW3rdoxWGThw8f4v79+20fBhE5Lo7jcWNky3ysNhtH3UGCesSFzYIgsCbrqZIMNrAmOCoznGZ6CMsFzCARUV29Xu9aY2TDajjXs0cutEs2ZWgm2XJs1gRHZTQxhFVGkiSt3BAqg+TrhEwiqk9lQ4rUUS6tVrUhkJvGlnZpniIZmizLWlkgtejY0jRtpL11507QLE1TbSnUKIpauyH29/cbWeZPROapOSo6O1wqG+JCo12GC4GcC3OOZgmCoHK7YvIzN/UAeWtK1+XlZaPvF4ahNeOsdahl/gyQiNyn5qg01eFSHURbszCzuHK8Lsw5mqdq9sjUZ24qawRYFBwtLy8v/B31YFFXI3FTVIB07949J3pTRHRFVfZt1Wmqg+havWHz8frUPlXJHrmUNRoMBjN/Zm8Jm0I9WNSlSDzLMsRxbLynI6XEs2fP+Cw2Ioeoyt6FOk3XENGs6Qy2LCevq8q1dPnRTZNcKMvK0tLSzJ+xFZ0hv8S1jiAI0Ov1GglY1LPYtre3nSqgRNS8ssGIriGiWdMZ6k5zaPsp7kXMOkbXH900z7zsjM38vBoaTC5xVZIksbqHk2UZDg8PcffuXQZIRDRTkWBETWWwPegA2n+KexFNHWOdNirLMq1DY/OyMzazpiSZXDKo82JHUWT9RO40TfH8+XPcvXvX+mMlInupqQx1GnQVWEkp+ezMHJPno2jgO61dDILAuY61iakr1gRHulJv0wqbixe7rjRN8eLFC9y6dYsBEhEZMWuYKE3TcV2sAishhBMZqKYIIa5Nds6fs6bev6l20fTnMjF1xZrgaNpqtbInVGWIeANeSZIER0dHuHXrljVbshOROxZl3GcNE4VhOHWV02RjbKrRtDlLNW/eUVPbsTR9bkzueG3qs1gTHE0rLFmWlQp01DNibB93blKSJDg+PsbGxgYDJCIqZVZmQVeDlB9y02kyK2MTdWymA5R5fz/frjYRKPX7fWN/21RCxJooYlpk6eoMftt6LXEc4+TkBDdv3jRaSInIfyrzoaOOU8GXCmRsqjdNaiJ4mxcwqPOuY4Vf29esaEKk7Oe0JvJoe0a7roBGZ8WhUxzHePXqFW7cuMEAiahhttQH846j6DGquSomGvdFDZgt59EF07J+k+dPx7wjV65J2aE9a4KjtncVLTuEN4sqbDZmvIbDIV69eoW1tbXWg1GiLrFliKfoUEtbwjC0/hhdZqKdDYLAiQBpaWmp1Oe3pgWvetF0pne7sKJtOBzi/Pwcq6urDJCIGtZ2IzKv06YCk7aPcVoApI6pC3W0SVEUGQkw2y4zRTkZHKmhnjIHr25kVy6MLQaDAc7Pz7GyssIAiahBVTNIReu4unWhDVMC8vNhlLaPySc6z+Xkdg06mbjmZaaULPxEQoifF0IcCiF+L/faW0KIbwgh/nD0/1uj14UQ4qeFEB8IIX5HCPHdZQ9+OBwW/l2bh7BsNxgMcHFxgeXlZQZIRJaqOmG2asNiS306GaTZcEzTuDjMpzP7ZjJobTsgLlLi/hGA75947ScBfFNK+XkA3xx9DwA/AODzo/++CuBn9BwmmTAYDHB5eckAichCKiteZmXT5KovFxtv4Or4VUBk62doetNGG5kMWtsOiBe+u5Ty3wH4dOLlLwH42ujrrwH4odzrvyCv/DqATSHEjqZjJQMYIBHZK78nTpmVZm03LHWp7IatgRHQ7A7TNmri2rQZfFa9g7allM9GX38CYHv09S6Ax7nf2x+99gYhxFeFEN8SQnwrf5Jdv6lNMF0I8wFSr9ezYt4BkQkulWsde+G4Xp/a/DBZW4+ry8re3/M2Rq790C0ppRRClK5xpJQ/C+BnASCKovG/5y7OV+napnskg8EAWZZheXkZWZZhOBwiCAJrliATddmi+9CW+1Q1TrqOx5bPZYssy6wJypo4DtPXf14wVfXTHajhstH/D0evPwFwP/d7e6PXCrPlwrdp8hw0dU7UKra1tTX0+30+p468w8bWPJ5jc3hu5ys7DyxJkpk/q9rqfh3Al0dffxnAr+Re/7HRqrUvAjjJDb9RQW3dAEIIDIdDnJ2dYX19/dqyR/VQXyKiWbi1ilkMjubTOdqxcFhNCPHPAPxJAHeEEPsA/iaA/wvALwkhvgLgIYAfGf36vwHwgwA+AHAO4Me1HCU1Qo3vDwYDSClx8+ZNnJ6eIo5jZFnW6cmHRLSYzQ98raONqQ5Uns5RFmFDlB9FkVxfX2/7MGhEDaVFUYSNjQ0cHx8jTVMEQWDVmDcRURPUlgq2SpIEUVR7CrHzkiQp/dy/k5OTb0sp3598na0cvSEIAgRBgDiOcXx8jM3NzfEy/8k5SBxqo7bNmzdApIPt9RyzWld0PhCZwRFNpR7Em6YpXr9+jc3NzfG+HqqiUNkkojaxYSDTbBhhmcfmrFaTdJ4Htmw0lZRy/CiBOI7x4sUL3L59G71eD1JKJElSaPIbe/VkGhsGe/i6RxqHrLqHwRFNpQKjLMuQJAmyLMPx8TFu376NpaUlRFEEIQSSJBn/Nw0rFaLu8HVvNNs/Ezuh1cxrnxgc0VxBECCKImRZhsvLS7x48QJ3794dZ4/CMEQURW8UMt6sROQCH56Rxk5oNfPmkjE4okLUvI40TfHs2TPs7Oyg3+/P7FGpm3VeVkkHHyo2ImoP64/uMrFDNnVMfp+PNE1xeHiI7e3tmROyVVA0Laukk87VCUTUPWqKwDwmO3hq8QvZhcERLaRWpak5RmEYIssyHBwcYGdnZ2pq0nRQRKTEcdz2IZDnoigyVs7U1ilVcPqCOQyOaKHJ7Iz6OssyPHr0CA8ePECaphgOhwB4w1Kz+LBqaoKN5YxDgvWsrq7O/BmDI6olDEM8fPgQDx48QBiG3KmVyCA+u4zyVMCmslrMopZzcHAw82cMjqiUaYFPmqZ4+PAh3nnnnfEqNlM4AZu6zNdnlzXNt3pEBUk2ZrdsNm84k8ER1dbv9yGlxIcffoj33nvPeOYoSRKvKraibH+EAZErfF3IoaY2UDFra2szf8bgiGpTQ2lCCDx+/Bj379839l5hGKLX63lZsS3Cx2RQnq+7UTfBlw7WZDDU7/e1/N2urKBjcERaqWEz9X+VKVL/VwESJ2Y3R0rpdGbJt2GOJvi6G3UT8sv3XQ6UdAVDk+qsoHPJ2dnZzJ/5/+lJuyLDZh999BHeffddVt4NkVIiyzLEcexkj8/XYQ6yX5F9jkxK09TJe9Z3DI6osHwmaNGqtCiKsL+/j729PeO7ZNNVT6/X66HX63nX42PZcQcb+vLUcyzJLrwiVFjZidZSSnz00Ud45513xhtI6sbhGP9xawh32NrQd2UODeljXykmJ5RpsB4+fIi3337bSIDE4Rj/MXNUDs/Xm7oyh4b0YWkho6IogpQSjx49wt7eHoIgYOVNpTBzVA7P15t0Zo5Yf3UDgyMyTgVIT548wb179xggOYbXyi28Xm/SmTli8NkNDI6oEVEUIcsyPH36FDs7OwyQSmj7PLExcAuvl1lt34+0mI55qAyOqDEqQHr27Bm2t7fH39N8bTd2Lu+fRKRb2/cjLaZj7yoGR9QoFRA9f/4cd+7cYQapgLbPD3fmpq5p+57rIp2rjnXsXcXgiBqnAqQXL17g9u3b6PV6zE7MMaunyi0MiMxgh6B5OtsAZo7IWSpj9Omnn2JzcxNRFDm9jX/T1I7YRKQfO2vN0xmQ6nj+JoMjao0KiI6Pj7G5uYl+v88ddgsSQji/bwsDYTKlbtnivKLm2RaQul27kvPCMEQcxzg5OcHNmzexvLwMIQSklGw85/Dh/PjwGcg+3A3bTbYFpAyOqHVhGGI4HOL169dYW1tDFEV8LMgCPmSOuvBUeZbh5gkhvC9XTWqqE2PbveJ27UpeUA19HMfjAEllkLIss+6mISqqbAaDZb0+HzoONmkqC2db2WcJIiuoymwwGOD8/Byrq6vjVW1dGX7hUAAx41FfV+qLpjT1/ErbAlq7joY6T2WQzs7OsLq6ipWVlfHKLFZ45BouCW8O6wfSicERWUOlw6WUGAwGuLi4wPLyMpaXlzuxdN22nhPV53uZtYkKjjjn6DO2rQBziV3Tw4nwWZA0GAwghMDq6iqEELi8vGz70IxSDalPQVKWZV59HiLqBgZHZB3V65NS4uLiAkKI8fDacDgcT9RWjS4bYHt1vQfPctkMzjOajsO61TE4Iiup7FGWZbi8vISUcryCbTgcQkqJNE0bmyzYhMmGVH0+l/lybchuXKE2nQ91SFsYHJE2anxb180ohBj/LZUxWl1dfWOTN18b4FmVva8VHjOAVEeZeiDLsk7MTfL985nE4Ii0Uc9LM9Fwq4xRlmVYW1tDEAQYDofa38cmkxVbmqbjCetV2dwo2HhMZB8dnQNb7wHd2NmojsERaSOEMLYFvFrKr7ILa2trSNMUSZIgyzIvMymT1I7S+c+q5loUrQRtbhRsPS6yi44Gv0hZ8zVDS8UwrCStTDVw+QpxOBzi9PQUN2/eRL/f97Z3NLkMV53byXNcJpM0+W+51Jxc01QQ7Wu94qIsyxrfloBXnyopU1B1FeowDMcVVpIkeP36NTY2NtDv97X8fdsUqZzrTkRlA0BdoBrXaR2JWfWT7iBMLSKh8oIgaLyuYs1IlZQpqDoK9eQO2VmWYTgc4uTkBJubm9Y90VmHopUzh6Oo6xYFHqpxnXavBEEw9d8mSaL1GCeHxItgMPWZpus564KjOI7bPoRWTfZsdN+gupQpqDoK9WTFppbwx3GMo6Mj3Lp1y8sAyTazKusmK3E2GDSpSOAxqx6a9W9tmG9kwzG4rmp9YV1w1PWNvPIbIAIYN/i2Bkk2SJIER0dHuH37NqSUPFcGzaqsm6zEq74X51dRGczIumnyPq9aX1gXHPk6f6SsyWiXWZH5kiTBJ598gs997nPjLQWI8ji/qhmz5vaQHqzb5tN1n7O2sBSDofLSNMWzZ8+wu7s7HnJjtoCoWT7tWm8jtg3zqSco1MXgiLwRRRGyLMP+/j729vaurXBIkqTz89mImjC5eIKoSf1+X888Vw3HQmSNfID09ttvI0kSJEkCKSV6vV7bh0fkvVmrwlzFYaxuYnBE3gnDEGEY4uHDh3jnnXeuBUas6IioDFeGsZgZvxLHMYfViKbJT8j76KOP8N57772xClBhsOQHNgxkiitli5nxK71ej8NqRLPkh9LUHKQkSdDr9caVXZIkzvQKaT42DESkE4Mj8lJ+KE1KiY8++gjvvvvu+OdxHC8MjLjdP1F5aZp6tUqUgXc3MTgiL+XnGKlA6PHjx7h//z56vd7453EcI47jqcNrVbb7J+q6/DMQqTmuDP+5giWYvBZF0XgoLY7jcYCkvleB0mQWiXORiPynht99wAyXXgyOqBPyFceHH36I9957b+5u7FEUQUqJOI6NDq2xt0dEZB8GR9QJajK2oiZpz1rVoHbX7vV6RofW2Nsjak8URQtXNpnswGRZ5tX8LJ8wOCLv5VelqaE0KSX29/exu7s7NTNkOigivZiBI1PyHRjd5Sy/i39VnAJgBoMj8t6sVWlSSjx8+BBvv/02siwbbx7GFWruYQaOmmBjOfNlzpRtGBxRp0ym0MMwxOPHj7G7u4soipCmKTNG5BQ2jt2WX3mb/z/Vw+CIOmVaFilNUzx9+hQ7OzsIgoATsMkpPj3HrElpmnoVWKogycbslosYHFHnqYzR06dPce/evbmr2OpSDVnVeQIc8iPSIwxDLwNLdsD0YHBEnZdl2Xjp/vPnz3Hnzh1jjxVRf7fq3+eQHwFXZdanrEcbfMkcTQZDdTJHXD33GQZH1Dkqa6MqFRVwhGGINE3x8uVLbG5uOt+rVBWd7Q2AC8domyAInC+fbQvD0Iuyp3MYTcfqOV/wLFBnqEpQZYlmVQJJkuDly5e4ffv2eNm/q7Iss76HXKSht/n4yV1tD60xU2MvBkfUGflKKEmSuUNUUkqcnJxgfX0dQRAYe8yAyUY/CAJEUVRoozvbsQFxhw/ZmKYwU2MvXhXqjHwwVKRCGg6HOD4+xsbGBvr9vpEMDBv9YtiAlNdWgGLLkJ+UkvcXVcYahzqpyMRmIQTiOMbp6Slu3LiBfr+vvVfMCdbFsJErr+vnTAhhLEhjZsx/DI6I5gjDEMPhEK9evTIWINnMls/JILK8rmfbpJTGym/XA88u6PbdQ7SAEAJhGCKOY7x+/Rpra2utBEhtBSlsBNzV9WsnhDAWINowbEhmmdnMhcgjqpIdDocAgNXVVQBXc5KaymhIKVupkLuefXAZs23m2JJRpdnqPgqKwRFRAdMCJCEEBoOBNRNQTWgrKCNqS5ZlCzsFDDztVzeAZXBEVFA+QBJCYGVlBQAwGAy8rSyZOSKipuh88HfdTh2DI6IS1BykwWAAKeV4iM10BmlWkFKkl0tE5AKbstQMjogqmMwgSSnH3zcZrNhUmRAR1WHTNikMjogqmJVBUnOSmuLinCBVAbp23GQ/tXy/TgeFmdj22FQnMDgiqiE/SXttbW08SVv9zPTNblNlQmQDriRzl02BKYMjohrULrxqSG1tbQ1ZliFJkkYCFxeDIxePuY6i2T0Xs4C2aXpYuyvSNG2ks2fTHEo7joLIYUEQjIfYXr16hZs3b6LX6yFNUyRJ0vbhUcvSNC30ewyM9OB51K+p7UpsunYMjog0SpJkHCBFUYQgCAo3jkrZ3ye7+brNg434sFkzmgpaGBwReSgIgvEcpOPjY2xubiIMw/EwW5m/Q/6wqcL3zWRHgsNqpAtLEZEmqhGUUiJNUxwfH+P27dtYXl4GwOGVrmIm0BwGQsXEcdz2ITiHJYtIsyiKxvseHR0d4a233sLS0lLh4ZUkSbybq9Tl4Q4Oq5nT1TJVVhRx7VVZDI6IDAjDcLys/9NPP8WdO3fGS4ynBT7516Io8q4y43AH6ab2NKLFmI0uj7UVkSFRFEEIgTRN8cknn2BnZ+daFmEyIJqGQzJE0wkhvOtEmMJhtfIYHBGNxHGsfThLBUNCCBweHmJrawv9fh9JkhSq2KcNyUgpvRt2U9I0ZTaAGtGlstbr9do+BOcwOCIa6fV6xnqiSZJgOBzi4OAAOzs76PV613pzZYKdyR6zT/N51HAkkSkqKGJZo3kYHBHlmMrIqHlEw+EQT58+xb1799Dv96/9vOjxqdVwStn5PF3qMZM/dN2bDIqoCAZHRDkm5zCooTQVIO3t7S2spCcbBDWPKT/cNhksLVKlceDcJ2ob5xd1W5ZljdZDDI6ok4pOUCwbeMyjKvderwcpJfb39xcGSEUahMlgyQQuRyfXqM1X5w05q0ysaT7PE2yKekxTY++36BeEEPeFEL8mhPh9IcR/FkL8xOj1t4QQ3xBC/OHo/7dGrwshxE8LIT4QQvyOEOK7yxyQesI5vfl0aa440KfoBEWTgUc+QHLJrGCxycqfDQ0tEgTB+BE+s6hMrGlcWdeOOvVEkcxRAuCvSim/AOCLAP6SEOILAH4SwDellJ8H8M3R9wDwAwA+P/rvqwB+pswB5edhdN3kTcsVB+7K36T5r6WUePz4Me7fv9/GYVUyK1hssvJnQ0NVszFl5tzZ0CFlR6C6OvXEwuBISvlMSvkfR1+/AvAHAHYBfAnA10a/9jUAPzT6+ksAfkFe+XUAm0KIncpHSJ2UZdk4HZ7/2lX5m3TaDasCJBsq4zJYcVNbZmVjFs1NKTPnrkyH1NS9y45AO0rNORJCvAPgjwH4DQDbUspnox99AmB79PUugMe5f7Y/eo2oMPUQ1y758MMP8d5777V9GKXM64FLKRHHMSdzU6OanpuizAqkdM5bpOYUbn2EEOsA/gWAvyylPM3/TF7VkKVmtQkhviqE+JYQ4luuZwXILBcDpXzqvmh2pdfrjTNIrmRker3ezGMVQqDX63EyNzWm6RVNyrysURMLJki/Qi2OEKKHq8Don0gp/+Xo5QM1XDb6/+Ho9ScA8hMo9kavXSOl/Fkp5ftSyvdda/iIFqmzl8pHH32Ed99915m9WIrM38iyDHEcOz88SnazJWvETJH7iqxWEwB+DsAfSCn/bu5HXwfw5dHXXwbwK7nXf2y0au2LAE5yw29ElU1OpHRlM8OycwaiKML+/j52d3edCJCKzMsIggC9Xs+5DOAiLpQ/ap4tmaI4jgvPhWor61ZEGw8ZFoveUAjxvQD+PYDfBaC6fX8DV/OOfgnA2wAeAvgRKeWno2DqHwD4fgDnAH5cSvmtee8RRZFcX1+v8zmIrJem6RsPnp0XOAkhcO/ePTx79owZF4tJKZ0IYqm7VIA0qyOTZdn4kSq2UrGK7nvt5OTk21LK9ydfXxgcNYHBUbeooKDow1e7LAgCbG9v4+DggAESUUPm1U1xHC/MlqpJ2LrqtzRNEQQBg3ADZgVHfuW4yQmqwnAxMFIPeW1qe4Esy3BwcIDt7W2re3VEvljUaZt8aPQ0ahK2rmGqMAzZOWoYgyPSqokx6zbnGqmHvDa5gi7LMhweHuLu3bu1AyRd527RajpW5OSqIp22IvPsdK9SK/K3hsMh58FpwuCItGoiu2EyvdzUhMSyS/XTNMXz589x586dWhk31QOtW4GqYdFZhBAMkIga1u/3ndtIVhfdQSGDI3JKmqZGM0emgzt13FXeJ01TvHz5Erdu3ar1KJk62wzkLZpMXjazlqZp7YCKAZn/bF1RZYuuPoJLd7vA4Ihapwp0kYIdhmFjD4s0QVXsVRvxJElwfHyMzc1NL561l7/mYRjOXbJbpHyogIxBkn/UvdNEdppDU3ab1pHS3S4wOKJWqQnObexj0QaVbalTwcdxjOPjY9y8edPZACk/sX0yQJpVwZXpGfq2nxI1t3eQKptkrzAMjd/jrEGoVWrSYtlhGFuCKdW465jHU0Ycxzg9PcX6+nqtNHob51EFRflrX0QURXOHVGwoD+S+LMvGWUzqLgZH5KQ2ApJp1OTwMpPEywwjzhPHMc7OzrC2tlY5QCq6JYHOc60e8VAlBT5vnpP6HG2XCbLbovKhMlTMHnWbexvNkBNURsJU6tPlYRO1o7LqodYxHA4BAGtrawCAwWBQKugo+v7qepoa2qiyy/Tkv8k3atwTimYputOyzWUoyzKn60AXMDgiI1RGgjfwm9Q50VX5qgBpeXkZWZZhOBxCCKF1cmJTq/iEEIUDpVlBkM2NGrXPhzpJ7ZhN5jA4IiO4o2uzBoMBsizDysoKpJQYDodOPW4gX9EXDY4YBBHg/rPtqjwzjGXfPAZHZAxv4OYIIca7466urgK4mpM0q8K1OavnckNHzbMpOGrqvrL13vUJzzAZwxu4WSpAOj8/x+rqKpaXl2f+7qJValWzfjqyhZxQTV2iewic9GDrRTRFUw+W1UmtmovjGIPBACsrK+j3+1ODDZuG3CYDtTJBtWvXiPzWZIeQO4WbxeCIaAYXG14VIA2HQ1xcXGBpaQm9Xu+NDRQXBUZVK/mq/27yeFw899SOrmaoeY+YxTlHRFPYklWpw5VJ2nWOp6sNI11n8xy6eapMxlY4p9MsBkdEU6jdm10mhMBgMICUEmtra5BSXkvFu9qgEE1qM+BP07SVuoL3rlk8u0Qz2JZhqUINsb1+/Ro3btzA0tLS+GeLdr2uOqdBR7rfh3NPzWj7UUJ1ghROxrYXgyNyjg2PDWlDlWBFzUFK0xQXFxfXHjWy6BEe8yr9eceio7LX3WB0YfJqkiSduC8mg++2AwwGN35icETOabsy1KFKY62CnLLU+RoOhzg7O8ONGzfQ6/UWnsNFgdOsY7Hl2uSPz6chiDRNp2bnqj6vzhXqek77jD5/bmqHPzUGdYYPFWGVxlplgRbJsuyN1WnKcDjE6ekpNjY20Ov1Zv6NJEm0HEub8sfnQ5lRgiCYeu5d/Yxpmi4sb8Bn19Pmz9mFDGVX2F27EXmqagVf5N+pxjP/u1mWjYdd4jjGyckJNjc3Zz7lPgzDQgGSzWw/vqpc/FzzAqAgCApNaLb1c+c/l+0dBiqOV5KcxX0+8MYKNGWyIVENkHo9jmMcHR3h9u3bU/+uD6v1TGCZqyYMw5nlyfVh8vzncvlz2KzuprxxHJf+NwyOyFnspZUzWXEnSYKDgwNsb29Pbbhcq+irVIBlqSFLKk9HeSoy/NY01+4TF80aSi5qVoZ87ntWfreOSJLEyhuSCNCT4Xnx4gW2traczxTNm0OlSxRFzp8nl81r5KSUzmf2fF5tOBgMAFzNe2xalQCWwdECURRVijon5ed8ENkkTVMcHh5ia2uL2TgyLkkSxHGsvdPpwiKBRXzKQk0GQWqPNbWViO2sKUmXl5dtH4JRQRAgiiKvCj/5I01THBwcYGdnx/kGhtqXpuk4AJrsEEZRhF6vp6XTSfbq9/utZIl0saYWXF5e1va3mKUhKi/LMjx79gy7u7sM4qmWMAzHAVAbZcnlITaf2q3JLFGbQ2tlWRMc6cQsDXXdotUdszoPWZZhf38fe3t7le4f2yYrVx26mbXRIpmhu9y4PMTmc7vl0tCam6WHiOaatrojvzFkFEUzN4qUUlYOkKru4m1KFEWVAqQwDGsvH6br5pULTnIn2zA4IipI7Snkasp+8vES87Kr+QCpDBv3R6o6tyWKolLZh7JB2KIgsu62ATYFqQADIGpW3Qn/DI6oUarA2lZxF6EafpdT9mVIKfH48WPcv3+/7UOxXpIkpYKwNE0XBgtql/Oy94r6fQYj1FVpmtae8O9/DU9WUQWWFbe+AHHWLtm6PH78GA8ePGhkk0VXTJ7vshVx0fJf9NEaVf42kUvKZFIX3QNqYvg8DI6IWuJSI/bBBx/gvffea/swCmkiiJNScnPYBjEwpyodhVnUxPC576flncgYVgrFuDhMp0sT83x6vd54iM321TRN7ZSdzxaVvU9n/b6Lc9ma0MQ1JT8VyRJNw+DIcvlKgYHSbC5lYSa5FNjVWebvs16vV+r+nPb7WZZ1Yi4bUVMGg0HlbQN4JzqEvSc/FQ3sdAzjzFq+X5RaxTZvo0jTc6CqaGIIrOz9Ofn7dQKjuteV7MLrqcfS0lLljhyDIyJHVN2zJ0/t31M3QHry5Anu3bs3tUFXq/naDpDymRnfH1WhVlFyHpQfdNynVA+DIyJLqAdyzqOjkZ/c76gKKSWePXuG7e3tmQFS20OdXcy0mggC1eTztoNdm0kpZ967VTOpOu5Tqo7BEZEl1AM5iwRJNsiyDAcHB9je3m49ECJzhBCIoujaNVYBEyeQX0mS5Fowns/g2dBRoPIYHFmGvbNqfDpvKkhyQZZlODw8xN27d400ANOCRJuGjrr6gGsVMHEC+ZXJ+9WmYVyf6sYmsWRbJt/AsFAX53LPzKbGvoo0TfH8+XPcuXNHa6Mw2RtXysy9Mh24qGfUkf9cDYJV3cgsXzkMjizmcoNPxZnoZdZ9aGrZ1TJpmuLly5e4deuWtqzXvPNS9JxNfg4TK+lsyhJ0TZMrI11fQcYsXzk8W0QeCoLgjcqwzOoXtVpmlml/J0kSHB8fY2Njw5phwSiKrn0ONf/Dxu0GinC5cTahyZWRLmQJmywfvpdFBkdEjqi7tFc9yLSoeZnLWccSxzFOTk5w8+ZNawKkaZ/DxUmyUsrxf5Ovd3nIpMlraXuWsMnslu9bDTA4IrLIvMomCAJrKqN5maU4jnF6eor19XX0+/2pDTqVp7Ik+QBXNVAqS9jlIKnu3l0+mMyUmrQou+w6BkdkPd97KHmLJhrbNG9gXm89jmOcnZ1hZWXF+krU5eBtcvjUpvLRpLrX0HT5bLJ8NZkRdS37WkY37yRyStnhIJf51LgNBgO8evUK6+vrWF5ebvtwpsqyzPp5JLSYyqpVZfK+mwzcXA3EfVE0EPanJiayVJnK0Ldlt3Ec4+LiAmtra+j3+9ZlAdV8la4E374xXZZ03If5wM2mst9VDI6ILFGlQvQhm6ECjziOcX5+jqWlJfR6PasaCCEEhBDjSc02HRstput6zXuIchWzGmBV3qg9RSfVMzgiMqxKyt63sfzLy0ucn59jdXXVyiE2l+cddZnpYWgTf59lzQ6LOqAMjog00F3hFa2UXRl+E0IgjmMMBgOsrKyMV7HZYtpKMMWm46T52go8Ju9Dn+YO+mpRObF70wYimsmlHqgaThgOhwCApaUlABh/X1SapkayavOGOrIsG2+jwEbPDuqamGZiCIzDanZYNLzGO52s4PpGdm3MJXC1sR4MBhgMBlhdXS09B6npc5xlGeeJWGjW9VCvN12XBEHgdP1Fb2LmiKzAxqcclTVyMThSGSQpJVZWVq5llOYxlTWaR51jFSSRHYQQc7NHbVwrlg+/MDgia9hUuTSVtq+qyL4uKiNj03nNGw6HEEJgbW1t/P080z6HiYApf97y85BsPY9dtSh7VITO+3ze+7rcmekqXikywvWl6C40hPljnDUsafPnUBmjs7Oz8T5I80w2LEmSaP98qhGbDIhsPo9dpeOaCCEaqavaGu6j6pg5IiNM9ZCamgNic2M4LVsy7Xht/gx5KmNUNIOkTNu8sW4mqen5RWma1t7dmcpL03ScFWzq3Nv0bERajHckGWGqgXH9USI6eqmzersun5fhcIjXr1/jxo0b6PV6hf7NtM+rGjpXMpeT5VlKiTRNmWEwTM0jy2cJm+DyPdo1DI6IGqSjlzr5sFElv8uzi41rHMc4PT3FxsZG4QBpkmp8XMnETGaqVCbD1uN3tWxNmrevVR1Vg/Kyu7MvekA11WfnHUjkKV2V8axhtMmhgizLxv+5kNKP4xgnJyfY3NwsvM3/NC730G0+9mnDUC7tt5Wn6zznA5WqQW0+k1VEGIYMkAxjcGQxFn6qS2UhmhiO1FVe4zjG0dERbt++7d1jVKpSAa6Nur4PVL6M1jkP+Tl0iwIl9dzCMmwdara1nWNwZJl8QanTcyZqUhzHWstrkiR48eIF7t692/kASQVFtg61dZ2JwLDIppJF31cFRbbeR7ZmwXi3Wca3gEgVehsLP+lTdY7QPGma4vnz59ja2up0YFB2DlKSJK3db64M37pgWjATx7GWv2MTIYSV7V53axxqhCr0Nhb+RdTKIfV/al6apjg4OMD29nanA6QyoihCFEWtBElquToDJDN6vV6lAInKY21DNIMa168yvm+alNKLbFyRZetZluHg4AA7OzsMkEpQQVLTmtpuw+QcLNs6Q/l73USW1iVpmjZyfVjTEDlIBWyuBkgqKArDsFDAk2UZnj17hnv37nV68i9dMfl4nzae4beIygR2XZPzpxgckRds6+k1waax+rIVd9GgKC/LMjx58gR7e3vWBEiLhji6WC6bYDKDaFtgpLR1r9u0KWkYhqWvT9Ed9ycxOCIv2FqhFeFDj7Cpnq2UEvv7+9jb2zP+XkX0er25lW8YhgyQyGlhGDo7j2w4HFYehmRwRNQyHT1CGyaN1/kcSZIUrnyllHj8+DHu379f+f1mqbKf0KIH5jJAIt2a7lCpIXxbMkhF9fv9yllmBkdEHpg2adzmjQPz8g8BLePx48d48OCB1mNRq610nzd1bRgkkQ5tzEHq9XpWDbGZxuCIFprXq+eyUnvZ/IyuvCrzj5SHDx9qzyDVOZ4if5v0aWrVZhzHxt6nasDRxhykXq/nRJ2iQzc+JdUSRdHMXr0Ly0rnBXAM7tynhthsmaRNzWli1WYcx+j1eoiiSHvmL03Tqc+ro/bxipD35gVwLgR38+iYa+RDgKgmaTNA8tu0sm5q1aZ6r3wdoTvzZzJLOY1riz+yLGvtmBkcETkmn4LXsUGlLQGi2o28ivwqtqYCJB1ZBM5BKqfJYUkfh0BVECmldKJTFARB4cA3jmOtK+oYHBE5ZrKn2dRKNdOVaf6p5FWoAGl3d7eR3riOxtPHBtgEneW7y89/U3O0hBDWdIoWKZo96vV6WjtGDI6IHKQ7e1REE5VpnewRcFX5P336FPfu3TMWIJlcrWPDlgw20lm+m3q8iY2mDUG2OXRVxLTsUZIkxu8TBkdEDqrT8NcNQBQTmaS62SPgqrJ/+vQptre3jQSNJrNSRQLdLmc+uqSpJfMq+MiyrLXAvGxdEkWR8Q4hgyOijikSgBRNY+fZ9DDcLMtweHiIra2ta5WoOsYiD7u1dT+XLmc+uqTpFWxBENQKOIbDYeWg3cYhPgZHRPSGKpvMqayHLcNCaZri8PAQd+/eHafl1bCCEGJukKQaJmZoiIrp9/vWdI50YHDkCFbSpEOZIbV5q0RmVYJNzX8qKk1TvHjxArdv3772eVSQNK937kqGRl1PW4JSqq/prKWuIbUqGSDdAZWutpLBkQM4SbMYW4dBbKJjTg/Q3hPCq0iSBC9fvsStW7esTN+XMVnxp2k63sHepqCU6qkypFan/qs7pFaHzrpEbVFQNECad84YHDkgy7JGGiNXn7yscJfZ9tlWftTxJEmC4+NjbG5uOhsgqd59vkIPw1D7EmaazrayDXx2TFmWsf7D1X1e5n64uLiY+TOeTQc0FdELIVjJUi22Ze/yGdc4jnF8fIz19XWnMl+KWlXERrAdtpVt4LNj0l0m0jSdO1oxLVC04fyU7Sisra3N/BnvMiLSJgxDq3rYQRBcO544jvHq1StsbGyg3+9bdawucz3rXIRtZRvQf0zqOi563tu0uYu+Tf1gcEREWtnQg1SCIHij0k6SBK9evcLq6ir6/X5LR+aXrmSdbSrbis5jklKOh+hmXU81zWPy571ez8rzUxWDIyLSanIYuO1NC6f1gAeDAc7OzrC+vq41g2RbZqFJXfjsOqY46D5POqddFAly52UJp2WPXA2YGBwRkTGqJ9qmacGREAJxHOPs7AzLy8vaMki6dh93URcyR1WpMmF72VDB0bzjnLfiddpih7bv/6oWBkdCiGUhxH8QQvwnIcR/FkL8rdHr7wohfkMI8YEQ4p8LIfqj15dG338w+vk7hj8DEVlK7XukKty2s0iTBoMBLi4usLa2piWDpD6nTZ+R2pUvDyqo6FL5iKLIyflIRTJHAwDfJ6X8owC+C8D3CyG+CODvAPgpKeV3AjgC8JXR738FwNHo9Z8a/R4RdVS+QbCtUVAZpPPzcywtLWFpaanW3wvDEEEQjPceom6afDC0zSsMJzM7XZk/tsjCKyavvB592xv9JwF8H4BfHr3+NQA/NPr6S6PvMfr5nxY800SdpxqJadVBlUBCZ7p+MBhgMBhgdXUVy8vLtf/e5Co52zJm1DyVObUp+GiqTNqwQWnZnbgLhbNCiFAI8dsADgF8A8B/A3AspVTvtg9gd/T1LoDHADD6+QmA26WOioi8M6tRqJNR0pWuF0JgOBzi8vISy8vLtQMklUGaxCCpGyavvQ1z76axPaulWxzHhX+30FmRUqZSyu8CsAfgewD8kUpHliOE+KoQ4ltCiG/ZWGiIuqzpOQJVKugizz4rW7dcXl7i8vISKysrtYfY8lx5TltbuhAw5gMRG4eY62qrHS+6CGLa9gPzlKqRpJTHAH4NwB8HsCmEUNvM7gF4Mvr6CYD7ADD6+QaAl1P+1s9KKd+XUr7fpciVyAXqqfVNvVdVReqOKgHS+fk5VldXtQZIQLeDpHkNmO/nZDJD4+PEfSFEKxOvVV1VNEAqqshqtbtCiM3R1ysA/gyAP8BVkPTDo1/7MoBfGX399dH3GP38V6VPJYCoA4IgMJ5uV88KM1k9VA1GBoPBtQCJVVj9zIDvAVDpOS2eBcpND9Gp5wyaemhukTBqB8DXhBAhroKpX5JS/mshxO8D+EUhxP8B4LcA/Nzo938OwD8WQnwA4FMAP6r9qInIOFMVXX5iahMNRNW/PxgMAACrq6vIsgyDwaBT8zMmTZ7HyeXpXZYkSafLhtJkWcjvfm/i3AsbekRRFMn19fW2D4OIGqKe32Tb6p1JUkr0+33cvHkTp6enpSZ0kllpmlqxCgr4rDz7xHTwmyRJrQdA5wOjOsd4cnLybSnl+5OvM9QlosbZuKx5GiEEBoMBTk9PsbGxMXUHYGqebZkam8twVaY3rAyCoFZnQ22XYez4jPxVInJClmWNTbzOU3OaXGhUVCV+cnKCzc1NJ47ZR/k5T/MeYUH6mOy8BEFQOXOk5irO2jJjUpUgjMERUYe5EqCUZWJZcRzHePHiBW7fvl1rOICqCYJgfF19LLNtamt6zbR5bGUCmaLloMr9yuCIqOOKzNtI07SVDFNVJieTHx8f486dO0b+Ps03ufM46WFLsCmEKBTITMsYxXE8s2xMfr7hcLjwPRgcEdFCtkx8XaSJfVaSJMEnn3yC7e1tZ86L6yafVUb12bz5cpVhsDLlosjcQQZHRJayrYc8r0fXZFZpXsXZZLDy4sULbG1tzc1Sufg0chvZNPnaJYPBYOb9Yus5FUJUWvhQZgfsIr9n59khovEuui4Iw7CxAKnX6xVKi5uWpikODg6ws7Mzs6EJw5ABErVmaWkJURTNDZJoOgZHlnNpngfpZ8MQQr4Mzqpgi84V0KXf7zf2XvNkWYZnz55hZ2dn5rWals0qel+rXcS7hg25PkIILC0tcRuKkhgcWS7f4DBQoqZNbtTW6/XYcE3IsgxPnz7F7u5uqdUzi+7nLMvGy5W7Jt+Qs7xRGxgcOYTLh6lp08qcTz1QKSWSJKmdnZFS4smTJ9jb2ys80XXR/WzqmVGuKVre1MaiZJfJIXD1WB7bMTgios4SQmjbUFBKiYcPH+Ldd9/VcGRUVtMPPjVB1xzDy8vLyv9W97YdtgyBl+V2SSIiqkntAhzHce0MUhAEePjwIR48eNDJuUK2cDWLpGuO4fLycuV/G4ah9lGKfPZoaWnJiewRgyMi6jw1oVxX4/Thhx/iO77jO7T8LSrP5SxSk9mjLMsamdM1mT2yYaHJIm6WHvKWlBJpmo573VmWOdkDJHNMZWSqNKizjiUMQ3z88ce4f/++jkOjDtGZPVo0PBYEwcI5XUmSaK+DXRhqY3BEVlFzQNREVPWAUleoCb6q98ehFf1s2jto0bE8fvwY9+/fd6KnTP5Rw2N15hBFUdTJTqo7rY5H2tzYz5VNBV01OTzD1UZmmAqQqjQCi67x/v4+9vb2GCBRa+rOIVIB0rT2YzAYONGulB0+ZHDUgjYrSVbQ5AsTgaeanK2TlBL7+/t4++23ef8tYHL+S5qmcwPfJubeuBBEzDJrTt7S0pIT5ZrBkSPavknafn8iHXRnj9QjW3QPIUgp8ejRI+zu7jo1TNw0k3toTXuSe1PvrbgQRMyTnzKQ58Lqs9XV1VIBEu/SlrR9k7T9/kQ6uJI9Aq4CJLWTNgOk2XRmcMoOk5rMHvnQIZ2XPXJBmblXvEOJaK58dqYLj3KoEhyVeVbakydPsLOz80Zgpybzd53ODE7ZBR0ms0fskM6nM3gcDAZTg+KVlZXCf4PBERHNlKbptUacz1abrsiz0pQsy/DJJ59ga2trfG7VFhZ8RBB1lRrS1mFpaQlxHNf6ewyOPGLL8mbyx7RhK1O96zRNa1eObd4DZQKbNE1xeHiIt956C1LK8SpHcp8L2T8VjNtGZ3atyETxefUNgyMisoKOitGlrRNUgLS1tVXps+sIJrtq1rJ0HVwIctV+clXNO3culcl5E8kZHBGRFYIg6NS8DJUxevHiBe7cuVM6I9elc6Vb18paXZMBz7zAPP8ctbrvU9eiv7e5uTnzZwyOPOJSr9lXLvWaquraTrkm5Ic1pJQ4Pj7G5uZmqccqsIGnJqitLfJ127znEFZduZZlmfYhyUV7W817/hyDIw/41iC7/Hl8b6yyLDO6DF3Xtbe9DE1Ovo7jGMfHx7h58yZ6vZ71x0/FuXAt5x2juudN1m2qs6B7PqPa2bsKBkce8K0n73uA4TLT+/Pomgti+z0xbV5KHMc4OTnB6uoqoijyZl6HrZo4hyY2FNVNSjn3XIRhaLxOFkIYW+hRdQ4YgyMPsKIkX+jagNHVIebhcIjXr19jfX19bmNhckJxVzRx/upOfG6CEKJUp6dIsGfjSriyGBx5oG5jsqjnQN3VdK/X5Z2jdZwrIQSSJMHZ2Rlu3Lgxc/5GkWEOF7IWRZiqm1wua20q0l4wOCIr1F2SycCIZvGlgTVNZyZHBUjn5+dYWVmZGiDNC4zUIzNczi7lj9vVz+CrIsNsZRYW2Mr+DRk6Ri3vbfr9OM+HpgnDsNVeoOkJ4LpIKbUPn6g9WNQjD4o83DPf2bF9OGeefD3owvU3wZWyb4Iqx21+fgZHFjMdKOUDI1WhMkiiSW01si5lrEydo6IBkrqXTQRpbehqUKC4nPXTRZ2Dtspzt0ughaYFJ6YaCWaMyBRdFbvuRtKlgEsZDAa4uLjAysoK+v3+1HOrXvMtqHDxeulie5BrMqNsw0R2v+4kzzQZuNQNlLrey6HP6JrHpruhd3n+1GAwwPn5OZaWlqbOQfItKFK62nkrez3V5olNPtetieX9pgOkeefLzzvKM75WfNQO08/ksiUjORkIlV2ybJvBYIDBYIDV1dXKuxC7Jl+O2AGbTZ2nuh3cMue4yXupbsdmVhA073y5W1OQVWxoDHVxNbtQVBPBiw3lYdoxuBAcqd2CJxsqIQSGw+F4iK0rAZIpPt3nKsNSJ9Oi5qwtOi9tPfC4Tp0y676fd744IZu0c32VhQ0Nu0kuX5syXLmOaZpey2qpr2cdv5qUvbq6CillrQd9ukTn9fQpMNJJnWN1fqbVFUEQNPJIkWnHVdWsOm9eB6MbtaQHXNpUy5VGaRbXj98WtjVAtt5D0xqZRWVwMBjg7OwMa2trCMPQyyEnk+XH9SFWU1RWeV52OR+8NznHyQTOOfKASzeyi8GFj41Lm2w7nzZPxq46zDkcDnF6eoqbN296seneJJP1iIt11CJxHF/7fx2LyqT6Wdsryuqa12Fyp8XtONM3c9d3yvaxsjRpURam7Q3cJgkhKj+A0gZqx+tJSZLg9evX2NjYMPbgzrZ0+Z6sMq9Hle+mAhYdQZjN7Km9qFW2rDAis5IkaWSZvU2BEfDZvB5XBUEw9ZwKIRDHMU5OTrC5ueldgDSNrRlAnarM56mzo3iVznEURY3eU1mWaRsaLxLY2VWDUaO6ni3qoiLPRSqiyHwDmxoxl9P/RYYE4zjG0dERNjc3nc6QFWFb4G2Cur+aqp+rdI5VYN7UMQZBoO0+LnKP+F/KaKYms0VtbFJGb6pyvWf11hZVMLoasS6XGTWcVuRcJkmCo6Mj3L59m50eDzSdmanChWOcpsgxMzjquKYmqqqG1PderY9czrp0TZIk+OSTT/C5z32O160BTc+7sW3FpYuBUVEMjjqOS1oXy2/KZ1vlZJoNnzeKok5mj9Swd9kgRwiBFy9eYGtry6qhTR81PceLAW9z2CrSzJUwdEU948eGhyE2rWuf1yZ1hsbSNMWTJ0+wu7vLzo8hba3W6mJHoQ28azSQUiKO40Z62ZM3ho4bZdZKGCJbdHE4tu4E1CAIcHBwgJ2dHaRp6v3S66a1tTKwi/dCG9giaiCEQK/Xa6SXPdmb5I1CRLNkWYZHjx7hwYMHbR+KFbIsm5mR05mRkVIyw+M4BkeOmeyt2DAnhIjslWUZPvroI7z33nteT6AtYt7+QTo7mrrncjLr95k0TRtp9xgc5cRxbLwQ6p7fUzZbZfNjFIhIv16vh16vh/39fezt7TGj0RCdwZGOITw1/cN1YRg2MkrD4ChHVSI+4+o0om6SUuKjjz7Cd3zHd3jRSHaBzo0gXVxQkmVZa2WVrWTDuDkbEbUliiJ88MEH+M7v/E4GSC0om7mv8hiRRX/PJUEQaE9YXF5eFntvre/qiTiOjaWemxjzZ9qciGbp9Xr44IMP8N5777V9KJ1T9pFNrgUzLlheXi70e9ac+fPz87YPYazX6xlbBaZr2fxkry8fEHEFGxHN0+v18PjxY9y/fx9JkjCL1BCdzwcjs6wJjlZXV9s+BKcxICKisj766CO8++67bR+GUfOW70/SGSRyOb/brAmOqJzJcViuQCOisqSU+PDDD71e5l9m3o7O+S1czq9fU8v4AQZHAJpZwq/TtMJR9Cbs4vPBiGi6/DJ/Fx414lomxrbl/LZbtN1AU8v4AYeDI52P7HBtCb8qHFWyRWo5J4MkIlKklHj69Cl2dnasDJBUUMTpA34TQiCKIiuSFfbdBQWpR3YEQdDYc810U096r6puJeZrGp2IysuyDM+ePcP29rZ1W44wKNLDhY0gVdveNmeDI6XJ55r5hJtBEtGkLMvGGSTbAqSu4kaQmZHh1EXnlK1ji1j5EJFthBB49uwZPve5z7V9KF7iRpDlBEFgJHM4HA7nv6/2d2yJii6bWLVVdiOvWUwW0vzxuTjkSETtEULg4OAAd+/eRRRF7MhpVPZcuhbMuGJpaWnuz7056yq6bKIgpWmqJQhTPYK6c4+myQdErqVRyR1sNP12dHSEjY2NcT1F1eTvE10bAZNZvEIVRFGkNeAIw1D75GjefNQENph+i+MYL1++xK1bt9Dv99s+HKvN6iiUHUZrW5ZlTh2vKWxBDSsyBKezIKr3YnBETfB9xWOZ3ZV9pDJGp6enuHnzJnq9XqfPxzyzzkuRSdA2nVNXMlum7037z4DD1F5Ckxdw8qLqLIiM+KlJJoZsbWsofA8AFwnDEGma4vXr11hdXbVimbWN6tTjZbM1bdwjNt2XgPl705rgyMf0vNrQavKmmXZBdRW8eY2VbYWbaBqWUzsNh0OcnZ3hxo0b6Pf7vE4aCSFKNfRttJddy6IyOGrBtBvBdKFzbdyb3KG77LqQ0u8iIQSSJMHZ2RmWl5cRRRHrlBnKnpeyWRCdS9uL3r8m5sbazJpaqOuT/UwXOm76SKb43rHpUm+5iMFggLOzM6ytrS1cDt0Vk2XElTKja1uatpg8fmtaS9ceKFhXUwUy/z5divqpOQy6u0UIgTiOcXFxgdXVVQ6x4c36vMpcvDYCFdc7zZ2YkN21m2vyRjAVuHTtvFLzTFSuNpVbdirepJ5peXFxgaWlpc5nkHTcA2WH4jjvyOykbGuCo0UrIGy7KHWpeUemP5PLvQIiU3yqS9o0GAwwGAywsrLS+QCpLpdWRk4L5NQ95cs8tM496lhduLaDhlkr1ly5OYhMauo+4D1XjxACw+EQQggsLy8DuAqYqLyy5bCNJx8EQTA3Y6Xa1yzLWm9j63ImONJ1om2rCG07HhX9q0Yjf3xVCnyapnx8CdVioqJV5dq2+89VKiBaWVm59j3pZUN9OmueUv6e8iEz63ZoV0HZ/SSaZMtxqXM0eTxVjq/pz6R6NdNSu1XG6H1JEbtMCMHr4IDBYICLiwsOsRlkwzPu5nVUfOp0WBMcDYfDtg/BS2VvpHmFu0qBbzq1qt5v2rFWORYfbvIy2q54p7G5Q0PXDQYDnJ+fY3V1FVEUvZFBsLF8ucSleUmusyY4ajtV6CvXx33L0h3cda0isjVD07Xr4LLhcIhXr15hfX0dURQhSZJxUOTDcEvbulant8Was2zqGU3Tnm3WJWxUqAyWF9IhjmO8evUKN2/eRL/fH5crdoLtodpHk1zev9Ca4MgENXGsqxW+rVkAshcbL9JBCIHBYICTkxPcunULURQhTdPO1sU2auKRUi7XJ14HR0C3e8JMv9I803qNXb5fSC/1LLajoyPcunULYRg6nUmYxdXPFASB1me0TeNyfcLW02NdHk6kxRY9rqDrk2d5/9QThuF4J+2XL1/i7t27Xq5iczk74nLwkhfHsfYsGIMjz3W9gaPZoiiaWzm6XOnr4EvD0YY0TZEkCaSU44zRp59+ijt37nhXrnTVsV2fH1tHFEXaR0oYHDUgyzIkSdLIHKB8ilcI4V1FROa4OjxA9gnD8FrwreYcvXjxAnfv3p0aBLha/nQNTYVhyIC8IhPnjcFRA9TYbhNzgEyPIZO/2Gsl09I0xbNnz7CzswMpJeI4HgdFLH9kEyuDI5VpIaLmsNdKTRBC4NmzZ7h37x6iKBp36BY9fJz0kFI20r7GcWz8PUyyMjhqYhZ9lzSxnwW5j/ccNSVJEjx9+hR7e3tI09T5htQ1TUzxcD3YtTI4Ir0494gmcfIntUk9EubRo0d48ODBeNm/61z4DEII9Pv9tg/DegyOOoINIeUtKg9dzzTyfjFLZSmzLMPHH3+Md9991/lMA8DsaxtMLOMHGBx5Z1bPhfNJKI/L+Ofj/WJOkiTjYTQVEO3v72Nvb6/Nw9KiauYoTVM+0aCiXq9nZLFT4b8ohAiFEL8lhPjXo+/fFUL8hhDiAyHEPxdC9EevL42+/2D083e0H7Wl1OoLk73uRWPz83ounHtE1B7Oq7kSRdG1LFEURZBS4vHjx7h//35rx6VjSKxq5khtmEn2KHM1fgLAH+S+/zsAfkpK+Z0AjgB8ZfT6VwAcjV7/qdHvLXR+fl7iUOwkhECv1zPa61aVSpWKVs09YpBE1Lw6925XtBkgRVHkxJwhakah4EgIsQfgfwTwD0ffCwDfB+CXR7/yNQA/NPr6S6PvMfr5nxYFctSrq6sAgIuLi2JH3hGzbtY64/P5CdpZljFQIqqpTMBT5t5t4uGgtmk7QCICimeO/h6AvwZA3aW3ARxLKVXLvQ9gd/T1LoDHADD6+cno968RQnxVCPEtIcS38jc/J0Jep25WE71Nda67Pr+EqC5Tk4mFEHOHW7Is87LOVAES537ZbTgceln+gALBkRDizwE4lFJ+W+cbSyl/Vkr5vpTyfY61zhfHsZHKN1/x+lrJEvli2j0aBIG3AYSapO3r53Odapd8vT5FopI/AeDPCyE+BvCLuBpO+/sANoUQKge5B+DJ6OsnAO4DwOjnGwBeFj0gX090HSowMjUerjJ3PPdE1enO7k4OqfkcCE0jpcT+/j52d3c79bld4XNgBBQIjqSUf11KuSelfAfAjwL4VSnlXwTwawB+ePRrXwbwK6Ovvz76HqOf/6oskZJYWVkp+qudo3M8PF/xBkHAlRK0kCsTiZuaQzc5F0h3dnfRkJoJtk1IllLiyZMnuHfvHusoTWyfw2bLgqE6pe1/BfBXhBAf4GpO0c+NXv85ALdHr/8VAD9Z7xDJhDYqXvpM3UaoTuVRdXdsVzbpMzmHbjKT0wYppbEhcBsnJEsp8ezZM2xvb7PO0sD2c1jkiQ5lp4EMh8Pyx2HDPJMoiuT6+nrbh9EJKmPESdikJEliZaNYVpqmjZTrLMusb2CqsL0cBEGA7e1tHB4eWpFZcJHtZbeN9unk5OTbUsr3J1+39yyRET48Z63rzwWr0zBM63HZ3CCW0VS5DoLAyvJXd7jE9nKQZRkODw9x9+5d5+uwttgcGAH12qcy7UKRTJLdZ4oKs30cWacwDL2eCLhInYbBlkm9VVdHtlnOkyQZH7MN53CS7Q2fDmma4vnz57h79671wZwPkiSxqm1J03TmlIQy7UKRB+/6fzd5Lj+pmuy1KBCwMRNhkgrSygRJbQ4JqCE7G4OirknTFC9evMDt27eNB0hduy/z0jS1brGOugebGFa151NTJTYVXLpOpXnTNF14M/syh2JWYzJrp+cymaxZZb2JBoyBkV2SJMHLly9x69YtowsF8vdl1wIlW5/3FoZhI8Oq9n1yIk+oBjUMw4U9XF09YBWEqP8mf2barEyQydWRXZ+D1lVJkuD4+BgbGxvGAqT8fWmqnHEDXjsxOCInzKs8WLF8RgUh+XS4qnxNLgFXwjBE048DiqLIm8wblRPHMU5OTnDz5s1C80hmKVJOTZQzKSWEEMxKWojB0RxNNCYmuHjM00xmQWZ9LpsmDNpIDV3NG8LSWWZUyrvJ+8eWybm+3HsuieMYr169wtraWuUAqWgdorucdSEwaiODrYMdNYql1EV1bdmo6o24rmjFYcP1mVzFVOcatDHx2ESZmXX9bF7xVZcv955r1NLstbW1a98XVbcO4XWfbbI+c+VcMTiaw4ZGtwobJ9EVkWXZuEF1rQFVWRKVnalSAahKpI2eVRtlxrVrXISr954PVEC0uroKIQQGg0Fj7+1jWdZlMtvmyj3iXXCkGqkupCuLciFSzw/B2H6s00ze8FUqAPX52w7KTZcX9bdVMOwjF+45Hw2HQwghxs/obCpACoLA6/Jcl4vZI++CI3LT5M6oNt046kY2cUPnK422gyKlqYrLlR5kFS5U/r5SAdHy8jKklLi8vNRe1qYNfftcnutyZZ5RnnfBETNGb+L5qEfd2KaDI1vYdjwumncObbzmvhkMBpBSYmlpCWmaIo7jQufcVFCbH3bvosmOnwttUjevFFEJ+dVeJv42UZt8Xe15eXmJ8/NzrKysYHl5udC/0b1U39dz2wWsmck4VUGYqChcr3wYHHWPbdfchV58FUEQII5jDIdDrK6uotfrLawvilybKgGUyU1QyQzvhtWoW/IT8Ilspp7PZhuf750gCDAcDiGlxPLyMrIsmznEVnS4s0ygMzkJuctDa67hVSLjVGVgaljK58qd/NFmOXU9w1rXYDDA+fk51tbWZg6xFT1HRScXT8swsa5yB4MjcpqLlc20SrPOXIdZD3XtKlvPRZsZAxfvE52CIECSJLi8vJy5k3bRrF6VrBFQfbGQreW5KVmWtfJ4IAZHpI0qwEmSGHsPHx7SOK1yrdNw+rJCU9eDPV07F1UbvzL3gmvnxAQhBIbDIc7OzqYGSEXPke7fa+rvuMrkfK15+2AxOCJtmtivx2QgoII7k72UNE2nHn/dz+RDBao20qsbILl2LiaPt0gmUP3ctc9qAxUg3bhxA71er+3DWajr19hknT/v+jM4Im1UAW5id+U6kiSZmt3Kz40yFSBN9oBMZtlcoQKiRQ/H9dXk5y3SUw6CgBN7axgOhzg9PcXGxkYjARKHvu007x7i3UWlpWm6sFFvqtGvMhQThuHU7FY+uDPV8Ew2hDauXmpaPmPUtcCI2hPHMU5OTrC5ufnG8790yz8yh9zA4IhKmxVcTP5OU8dStsIpkqZtqpH2ORiQUo6D10UBbBiGXp8LMkuVtbLiOMbR0RFu3bplPECqku1rYyKy79QDihdhcOQYW4ZhbAkuAGZfbKWel+fLhPG64jhu7b1tqTcWqXqck89mnCbLMiRJ8kagniQJjo6OcPv2bevqEtuOxwdFh1EZHDlG9W5cqeyI6IqqlNsIkuZlRdoM2iZFUaS1bpNSIkmS8QaPURRNDdSTJMGLFy9w9+5dBiSeK9pRY3DkKNMpYJovy7Jr/5k2rcGoOpTgMh8+r20rpGw7Hp11mxACURQVfizI8+fPsbW1ZeVkd3aIm2VfCaDOUCluF296NX+giVVDSZK80WCo1S9d6uXa+vgN8keapjg4OMD29rZVAVKSJCz7DbPn6lPnqBR3lZ5ifpXatHkEPpl2ftQciy5ljxY1Dl05D2RWlmU4ODjAvXv3Wp8rp8r0rOFA36l5Ym1gcEROyq9u8qHiqDrvo8hEVBuoRwAUWblWlQvngfQxmXXOsgxPnjzB3t5eq3VLGIadDvrVPDET5u2ODTA4Io+pijOO46nBh00TUW2b96FbEARGV67ZdC2pGVWzzkVJKbG/v280QCoyZ5FBvxlLS0tzf87giGZSPf0kSRDH8dQlsDZTFWev15safPgekHRFHMe8lmREPkAyIQgC7p5tKQZHNJPq6UdRhF6v90bvKYoi9tjpmjYC5jK9+i7N0eqCJhZ0SCnx+PFj3L9/38jfD8Ow8uRv1r/mMDiiWlzvsbNy0SsMw8YDpCiKCgc8XZzI7jN1PZsocyYDpKp6vR7rMEMYHFFjVAVm03CcC8GdS5Wf6Xkgs5Sdl+HKRHZarEiZS9P02qKAqooGSE0G3y7UYS5icESNURWYixtYqspOVbBNYu+QqDx1r2ZZNn6YtI5n+KkAad7fUe/H7KS7rAuOLi4u2j4EKzTZGKoKhGZTmQZVwTaNvUOqostBtbpXTWzUWmQVm6nsJOvqZlgXHFGzq2/UjWbTbrBEVI8KihhUm6FWse3u7ja6D5J6RhyZZ91ZXllZafsQWtdkhdZWJoSaNW2uhc07i9t8bC6wLShy+VFBs0gp8eTJE9y7d09LwJJl2cIyz8CoOTzTRJYxMRQyba6Fjp3FpZSlj7fIxnc27XrOQG26Mvvz1HlUkKlj0vV+z5490/IstiAIrCnzxOCIyBqqV21br3+eJElKH6+aA1Kkp9y2JEnYaGF60KEmHdtE7cDeZICknsW2vb3NLLxH7CrZRB2mq1ddJDNTd0mzUieQmxV06Do2HaIosi4AaIPJQEj3Fh9tBG1ZluHw8BB3796tFCC50FHoGufuejVuzSWSbz6FnOeEABRanaNjSbMpNh+bbj7NwanK5S0+8tI0xfPnz3Hnzp3Sn6Xp7CQ3QV3MueBIjVszffnmxnc8J/7IV1y29yiLNvBFMlptauP4mgoIdJehNvb7ssW8c5mmKV6+fIlbt25pvba6r19+p3ib78k2ORccEXVBPtC1aZhpmqKNgIn9ZnSy/fjqmLeXWZWy5eoqVx330qIhsCRJcHR0hFu3bmmbP2iqDrBx3pgtOn1WpJRWNzqzTB6zi5+BiouiiL27DjFxP8+bO+X7JrD5jKCOIdswDBeeryRJcHx8jI2NDfT7/VrvB7AOaEOngyNXU8OTNwlvGn/Mahht7aUXachdCN5tOcY2hjmKTDq35fxUYSIjWOR+jOMYJycnWFtb0xIgmaoDOLQ2XaeDozAMnZwEyLlG/phsdFyrpBYdb5HsbNsN7+QxtplRtvGBuC43nouupenrHMcxzs7OtAVIJkwOrbk6oqKbe5EBUY6U0umVTWmaXgvQbWsYF1l0vGrfmXlUZdzW3If8MapAoMhxV+VambUxYCtKNfJCiGtf539uuuwNh0MAwNraGoQQGAwGxt5LB5fKpkmdzhz5qktRv6s9WsXFzGUV88qkapjaKreTWSNTk1SnNc5tcTkbVEZ+SG1a+WoqIB8Ohzg/P8fKygqWlpYaeU9b2NYeFT2ebtTMHTPZM3Wtp1qGqz3aLpkMPoA3A4S2skbqeNT/TR6HTfdg0eHONjN6us36HE19PpUxWl5evva972zqFACftYeL2kXrg6M0Tb3evt/EU5Yn/55thZO6xeZylx9Oy/O5QwEUCwhs6/H7QAVEKysrEELg8vKy5SMyTz0qyNT9VHYYPJ9JdDo40kWlkG3sBamLpDtQyg8RsKIjG9gacEw7LhMdF5d0+bOblA+QAHQmQDIpTdPS2zQsOibrgyPfh02aqoBsbZTatKjx86Vx7ELmUFWOuvh8rnxkMtNnoh4YDAYQQmB5eRlSys4MsZlgqo62PjjSxYVGTvcxmq7gF1UauhssWqwLgdCkNE2NfN4unUOXmRyyUUzUZZeXl5BSjjNIDJCqM9G+2x8xUCk2bWppe+Oy6IZyIaAuwuSydBPKrqLiIxCaoR5WamqVW9W6y3T5LjPntey5GQwGuLi46OQqNtuxRvHM5E1sMljqSnDhEtcCoWmEEKXKLctZM1TZmlwJW+a5X/Oua9Vy20R5XzT8XudYGCDZibVKQa48amTyJmbDQUXYNGHf9UyQTedSt2lZlDKZlXnX1dVrnv/sVQO1wWCA8/NzrK6uatlJ28a2Sh1THMctH0kx1pRGtYuorVytsF3PIlAzJrMBbXO53M7aHsBHZTOVLl/XWXR9puFwqO1RI0EQWBcgBUGAJEkQhqETAZI1rb3tOwX7MFxBVESRct5Exatzbksbu0HrrC/yT5afpszQFtlLBUg3btxAr9er/Hds7Myrx9AEQXBtcnuSJC0e1WzWnD3bLiQRzdZEz7Ts3KNJ+WDC5foly7KFO1WroS1bGxoqbjgc4vT0FBsbG7UDJNuoY8qXZZVRso27NQbRiG3p4y5oomda94GnLgdEeZM97WlUo1PkfCVJ0onnqpUxrQ5ps16J4xgnJyfY3Ny0flSlrnnl28QzAIteVz9qD+o0H/ZSsrHntIiNPdOuK3JN1NDGPHEcextAqUxc3rQ6JAzD1gOko6Mj3Lp1y/sAaVa5VdNZdA4Zq2u96NoyOCKywLzKj5mxKy4GkDZRk2CLBFBRFHmTeZtUZHWdKmtVO15SSi3lNUkSHB0d4fbt2150AquYtn2EjmBp0fn0s/QTecTGlSdtiKKIAVIN8+avTK4e6nJWUK2oqkMNCeu4b5MkwYsXL7C1tdXZAElRgVEje1sZfwcqLY7j1laeTEs5U7vqzr1xXb7h9n14QZeyjXKdib+T7+v6cFwURVp2xNZ536ZpisPDQ2xtbXmb0SuiyVV43T3LFuv1eq3tA1FmQzcbMbPgnzbvB1epRrnp+yEMQyOTaG3U9IOp0zTFwcEBdnZ2Oh0gNYVn2FI6dkntoq5mFrIs82YSbZIkb2Qvp2U2GAgvpvt+KDKXpsiEbx+08RmzLMOzZ89w7969xjuxXRtV8L8EE3VAEATo9XrXKuw0TZ0MmKIoKrSpIecgNU8I8UbAlaYp58Q1KMsyPHnyBLu7u40GSEEQdCpAYnBEtSRJ8saQh+2PgumKMAzfCJhcUXTeR9nMiHpQKukThmEn5sTZNLQrpcSTJ0+wt7dnJECaFQSFYej0tIsy3Ks1ySpRFL0x5OHSkOC04M6mSpD0UpNkGSRRWbbNfZNSYn9/H3t7e1r/pppL1ZUgaBYGRzTX5PyPWcMYNlUaZUwL7nSt3CH7qCFGV1cAJknCoK5FttUNUko8fvwYDx480PL3bHwmWxllnzG4sbEx82fungVqxOTwxqxhDNsqDaVoVsjV4I6Ka3p1kQlRFL0R1LW59QfZ4eHDh9oCJJeFYVhqXtTJycnMn7ldU5A1VHBh23yjolkhW4O7PNvO7Sy6dgfWzfXACJieOer1ep0fAiEGSIqueVHu1xZkBRVcuDTfKE8FdzZnkHq9nhMBks7dgem6WZkj11YkkhkqQGKwXB+DoxKSJHGicQK4B0xZKrizOYMkhHAm+HR1To8Ol5eXjb6fqysS6U1l58xM8+jRo1qr2HQ9u8x11txR5+fnbR/CQlEUOdM4Tc4N4uocomYsLy8DaD5IIveVnTMzTX4VW5GgefK9Jh/02lXWBEerq6ttH4LXXO3J5zeY8+G5TWXYPMRHVwaDwcyfLS8vM4NLpemYM6P2Qdrd3V0YIDEQms6a4IhomvwGc115LIGLTGUlbc92BkEwt5ff1cfZUPvUTto7OzuFOsZtPhPPxukqbGmIDMovsy6bCbJt07l5TGUlbc92Fpkkz+wRtSXLMnzyySfY2tpaeC+1uceRjdNVGBwRGZRfZm3zZG8ddPc6bc8aKYsaFJuyR5xo2z1pmuLw8HBugNRm1sjWDiCDo46xdQ8ams6lgEp3r7Ns1kiV66YDAJeuUZIkCxtBmwOoLs051EkFSHfu3JkarLeZNZp3/7TZXlkbHHE5oRnqqdpZljnTMycqQlX6LNezFVn2XySA0qFMHa+Oh3MOq0vTFC9fvsStW7ecCehV0NbGPW1tSUuShBmOBYpWLNOW8QdBYP18DqKypJRWDWMtYmMHsKl9k9I0LTyk4lJQZOM1VZIkwfHxMTY3N60OkPLnsK22ytoS1+v1tFw8nzNQi4JH9bldXcZPVFaTHSod2RXfh4nmzWVxad+4Mmy/pnEc4/j4GBsbG1aef1vabGuDI12yLHvjyfK+WBQ8Tt6kPp4Dory6Haoy98iibEaRRtLlTkuRc6WG77tU94RhaP3njeMYp6enuHHjhnUBUpvzn/LaPwLDwjD07sGMRXsm+Yq3zdUIRE2rUtbr7kw8qeo8CdsbVqB4797H+reI/Llp63ouet/hcIhXr15hbW0NS0tLDR3VYraUf++DIx9Vqfg5tOYmBrTVTJ63IucxCIJKj1uYpUoWq86QQpONii29e1vlz01bjf1kQDrtOIbDIc7OzrC6utpagDR5XLPOV9PDbe7MXKQxlyacUnW6MxldolZkqkZK53k0mQUp8lwrKWXnMjEuazuIVOVlVrkZDocQQmBlZQXA/EfimJBl2bVge9r5UiMfTT73ja0s0YR8o9omKSWzfTXkA6LJ86h6oa4+ZHNaQ+Hi56D6FgXL6ufz6jQVEBUJkHQH5+phu4s0na0s9E5CiI+FEL8rhPhtIcS3Rq+9JYT4hhDiD0f/vzV6XQghfloI8YEQ4neEEN9t8gPQdByOcR8Do3oWnT8VIOlO15vek6WNIIj1ST0mM8DT/vbkUvgiBoMBLi4usLKy0vgQW5GH4zbdYS3zbn9KSvldUsr3R9//JIBvSik/D+Cbo+8B4AcAfH7031cB/Iyug6Xi2hyOcb0itSFrROao1VOmHmTcxL2nHnjLYddu053lLhIgqSE639U5q18C8LXR118D8EO5139BXvl1AJtCiJ0a72O1NE1b3ZF3ViDCrAPRbJPDUUXm+RTl45xAGzsMLu2EbirbN+vv1nm/pjJItgdYRUu8BPBvhRDfFkJ8dfTatpTy2ejrTwBsj77eBfA492/3R6/N1fQkMJ3aHOtX721TtsbGipRICcOw0j1rW2UeBEFn5xklSdLZz55n6hwsCpB0va9t91Re0S7O90opnwghtgB8QwjxX/I/lFJKIUSpTzkKsr46+trZBrXtDI0qpL5UFGpVQtvnleyVpmnj5aPqUEKSJIUzSSz7xRXddoGqUwmL1dXVa9/r0tTwXNX6olDpklI+Gf3/EMC/AvA9AA7UcNno/4ejX38C4H7un++NXpv8mz8rpXxfSvl+EARWP+fFBb4ER4Dbmaeu7QbcBiFEK0MqVe6xIAgKPdJEBUYul/0m8Tw1YzAY4Pz83Ng+SKbbrTRNK7/HwhImhFgTQtxQXwP4swB+D8DXAXx59GtfBvAro6+/DuDHRqvWvgjgJDf8Rri6YDYNg9nE1aXVSpeHOpqisgZFH1raJnWsiwIklT23vexMe4g1Na/J9iMfINn2qJFF6oxKFcn3bgP4V6ObNgLwT6WU/68Q4jcB/JIQ4isAHgL4kdHv/xsAPwjgAwDnAH680pF5LAiCcURre2VI5fB6NkPt+B7HsfVZ53zQM2+YzYWyM62xUQ01sznNUdnTpoZgB4MBpJS4ceMGXr16heFw2Mj71lWnTAobhgCiKJLr6+ttH0ajuMstUX11hqLUDuRNzvEpet9nWYYsy2YGUmmaWpNpUm3I5E7HZJaJNkRdy1l/t9/v4+bNmzg9PXUmQFrk5OTk27ktisZYiltiQ6VG5LoqDbGUEkmStDKpd9p9r44nLwiCuUFbEATWPF4mP7GWgVFzTLQhajRj1h5aw+EQp6en2NjYsD5jWxdLMgH4rKdKVJft81LyD2G2oZMy66HQk8eWn6s4OUcpjuNWA6VFwVxVLswr89G8KR/D4RDHx8fY3Nz0ck8vhcFRh6kVMgCXxpI+quG2PUCySZHjmdzRWwiBJEkgpUQURRBCtB4k6TQcDr1ufF0WxzGOjo5w584db68RW0MDsizDcDicuUKlrSxNPhgC3N5fiuw2KxvSVaY2uVVBkQqu1Pc+6PV63nwWHyVJghcvXuDu3bte3uvWtIyvXr1q+xC0Ufs2zYqo23oukgqGJoMkIjJLLYE2/SQAn4IJnz6Lr5IkweHhIba2trwLkKwJjnxbrbboxlaPMBgOh40HKswYUVdMW1GTZdnCuSy6Jzur+qDpp503hQ/B7a40TXF4eIidnR2v2hVrPklXewn9fn9coNReEkSkR6/XeyMQUpndeUGSTavBXDBvhRP5L01TPH36FLu7u94ESM58CiklLi8vvdlbYZqlpaVaQaLNE2CJ2iCEmLnkeNFji6o+oLYsKaXTD94GcG0VHXVTlmV48uQJdnd3nSkHa2trM3/mTHAkhMDy8rKW7csHgwEuLy+96eWooMi3MV+iLhBCvDHclqapUx1BlS3wpU6larIsw/7+Ph48eOBEgHR+fj7zZ1YER5M31OXlpdH3W1pawvLycuMXL0kSI9kdBkWkg0uNse/CMHTuOVZld+zmHkZ+klLi4cOHTgRI84J5K4IjIQTOzs7G3y8vL2v5u3EcWzXUFEURAxmyVr/fZ4BEWql9mKbxfYflLssHSK6yIjgC5o/9VdXr9ToZjLg+f4Ha41q2oq4sy2buR1YFg8vrmtx3Sfe1NKFL26hIKfHxxx/jnXfeaftQKrEmOCJ98vMXpJSssIlmCIIAURQhTVMtWWYVXPKea566lsBVoGTTqIGiNuzsSoAEwNkAyargyMbCDHw2gdtFQojOZQOovK5nG8Mw1Jpl5j3XLlPPetNh3j5zaZp6OandxQDJmuDo7OxMa2HWOflZTeAm8lU+29j1QImoLWEYeh8g2T5JW7EmONI954iTn4mKyw8D+bqLcxMGgwGDy5b4MJ8njuPSq/5c8vDhQ7z99ttOfD5rgiMiak9+GIjzZapbWlrC0tISBoMBz2PD8sNVru7W3ev1kKap80HeLFJKPHr0CHt7e9bvpG3N0blYkAeDgZa9Ony9EchNnC/zpiRJSt2nS0tLPI8tUhOfXdTr9awIHEwFmFJK7O/v4969e7U+p+nH+7R/BUZev37d9iGUtrS0pGWvDvUcJyKyk1rR5mInjqgKFVyaCpCePn2KnZ2dytNfgiAwmiG0Jji6ceNG24fQKht6CkQ2mbeBYJW/VVev13M2G0FuKZupLGLeKIeUcuoCJpMZuCzL8Mknn2Bra6tWgGTq+NgiE5GVVLamqnxjoPa/sR0zUwRclVddw0Zq7tu8UY629l9K0xSHh4fY2tqy7h5lcEStkVKOswOqZ+/DihPSp06FafPjKaY1ei7s8EzN0bW7eNG5b0EQtDKCkaYpnj9/jjt37mgLkHQElgyOHORL71IIMa4A1E0xb4M026hVJfnVJaYnCVL7dFxfIcS1v6OGNdoM6Fwoty4cI5WXJAlevHiBW7duaZvHy+Cog/JZFmpPGIbjnXhVQOfzHiV0RWU7695/+XIihGg906Vzjpdu6riYVdZHnVNbsvVJkuD4+Bibm5vaAqRa/772EWjW9bRykcpJFZyunyuiNqiJ2U3df03t2dPr9aytU1TjzY199ciXKZuy9XEc4/j4GBsbG61vhWHHGck5Oztr+xBaNZlqn6dudG1Lj4HIRU1letSQW1MBko0YFOllU0A0KY5jnJ6e4ubNm63u1m/d2bH15mxKvsCargx92FGWSPE10FcTZXV9Pl/Pk242nCfddbINn6mI4XCI09NTrK+vtxYgWRccra6utn0I1mgysnd5R1kiANcmxftG51Pmi2yP4OM5VIoumrDhHOjutGZZZsXnKmI4HOL169dYXV1tJUCyLjgiKsqVm7yqJjJ5Pp1DteLRp880Td0Gs0h23udzWPTc2bDvju4OchRFMz+/jaMHg8EA5+fntQOkKp+LwRE5y8abuYwiDZAtn8+lxtKGRs2kWWVC5zXy+RyGYehUllz3sU7LQKr5p7bUN3kqQFpZWdGWQSryORkceaLOTsKucq2Sm7ToBm3isxXtmdpSaXaxnE+atV1EnWvk2nl17Xh1aHIOqm0GgwEuLi6wsrKC5eXl0v9+8n5hcNQBXawkfFFkDoktwZ8Nq4VY1udbdI18On+zngXmK9NZcpsDI2UwGODy8hLLy8uVAqS8Ip/V7rNBhdnQeFFzutQw5BUt56b260nT1Mtz70L9kR829HnYbxbbg5cmXF5eaguQFuleCbNQmqaVd1Z2oVIj/WzJKDWpTFnXdX7UvZllmfO7n9etK9Q5aMvkue9S3edyuZulanm6vLwEAKysrEBKicFgoPvQADA4soIQAmmaOj+HhprDXuR8uhpOdW/yHp0+b0NK2VhZ7Op595UQolaAJKUcb/1jIkBicGQBVThM3fxqfL6LqWiiKlRAlL83pZSdbqAnP7vOcxHHcac2AFZDhF3u5NQtPyogWl9fv/a9LmwtLWHiJlEVjguT7YhsEgTBeN6S6c6Ly8qckziOEUXR1H/TpSEy4LOMJMAAqQ4VEN24cePa9zowOPJYvsLp8g3oE2YB59OVgRBCzD3HcRwv/B26bt6QZJn6SQ1vukwI4fxnsEU+QJJSYjgcavm7bDE9VjYgUnMIkiQZp32TJLFmjxtipbrIvGBFSok4jgv/rXmP1ImiiNehJF0dNF86ejof2eTjCsoyBoMBTk9PsbGxgX6/X/jfzdtU0o9SRlqomzUMw3EF5NIEVBXQ5XsOZRpDV7hyPdow69yowEhXpofPImwPz/ubgiDofIA0HA5xcnKCjY2Nwtnjee0DgyN6Q77yMVERqYaqbE9+ERXQ5W8MDnvU50OlK4QYz7+zganlx9RNPswr1fH4m+FwiOPjY7z11luF6v557+n22SQnqbkapuZsmA7u0jRFkiTebgiYlySJN8NHtgRGwFUAzwCJdLKpfFehlvbXFccxXr58ibt379ZqXxgcUSvUjeziDR2G4XjOiY7AwdRuzjpMq1xsaNR1TbpsSxAE1+ZGqLl+RDZK09T4w6d1Zr+SJMHz589x9+7dynU0gyPygssZHJeG/gaDQakJjyYMh0Mv9sSZzHDmtw8gssXknl+uSJIEh4eH2NraqhQgufVpiWbwZeinaWruV1FLS0utZ/v6/X5rx1D2fJURBEGhQDm/mrRLXM8Wuiq/QKctVR+8m6YpDg8PsbOzU/ozMDjqkCzLrBgSIXvYNEnZBkmSzM3eVD1fcRzj/Py8zqGNRVHUemPVhrYzltSeOqtD0zTF06dP8fbbb5e6b7p3h3VYEARYWlpCmqbshRFNEUWRkfk/vV4Pq6uriOMYFxcXWv82Ec2XZRkePXqEBw8eFA6yGBx1UBiG7IVRY7Isc2ouTa/XMzYPrNfrYWVlxcjfLkM92ZyuDIfDRja7lVKWGhJ1eS6lbbIsw8cff4y33367UIDE4IiIFqozHFt0Lo1LmmpMyygzH2p5ednw0ZhjIuvd1Dy2shObOZfyM2UDy1l/49GjR3j77bcX/i6DIyJaaN42+2WoPaJc1+ak8Fl0zh+TUlo7P7HJrLeJc8DsUTW6VsxJKfHw4UO88847c3+PwRERNUbtEUV2E0JoC4hdZuIcMHtkh48//njuELewITUshHgO4AzAi7aPpePugNegbbwGduB1aB+vgR18vw4PpJR3J1+0IjgCACHEt6SU77d9HF3Ga9A+XgM78Dq0j9fADl29DhxWIyIiIsphcERERESUY1Nw9LNtHwDxGliA18AOvA7t4zWwQyevgzVzjoiIiIhsYFPmiIiIiKh1rQdHQojvF0L8VyHEB0KIn2z7eHwmhPh5IcShEOL3cq+9JYT4hhDiD0f/vzV6XQghfnp0XX5HCPHd7R25P4QQ94UQvyaE+H0hxH8WQvzE6HVeh4YIIZaFEP9BCPGfRtfgb41ef1cI8Rujc/3PhRD90etLo+8/GP38nVY/gEeEEKEQ4reEEP969D2vQcOEEB8LIX5XCPHbQohvjV7rfH3UanAkhAgB/N8AfgDAFwD8BSHEF9o8Js/9IwDfP/HaTwL4ppTy8wC+OfoeuLomnx/991UAP9PQMfouAfBXpZRfAPBFAH9pVOZ5HZozAPB9Uso/CuC7AHy/EOKLAP4OgJ+SUn4ngCMAXxn9/lcAHI1e/6nR75EePwHgD3Lf8xq0409JKb8rt2S/8/VR25mj7wHwgZTyQynlEMAvAvhSy8fkLSnlvwPw6cTLXwLwtdHXXwPwQ7nXf0Fe+XUAm0KInUYO1GNSymdSyv84+voVrhqGXfA6NGZ0Ll+Pvu2N/pMAvg/AL49en7wG6tr8MoA/LWx7doiDhBB7AP5HAP9w9L0Ar4EtOl8ftR0c7QJ4nPt+f/QaNWdbSvls9PUnALZHX/PaGDYaGvhjAH4DvA6NGg3n/DaAQwDfAPDfABxLKdWD3/LneXwNRj8/AXC70QP2098D8NcAqAeN3QavQRskgH8rhPi2EOKro9c6Xx/xIUc0JqWUQgguX2yAEGIdwL8A8JellKf5TjCvg3lSyhTAdwkhNgH8KwB/pN0j6hYhxJ8DcCil/LYQ4k+2fDhd971SyidCiC0A3xBC/Jf8D7taH7WdOXoC4H7u+73Ra9ScA5UWHf3/cPQ6r40hQogergKjfyKl/Jejl3kdWiClPAbwawD+OK6GCFSHMX+ex9dg9PMNAC+bPVLv/AkAf14I8TGuplN8H4C/D16Dxkkpn4z+f4irjsL3gPVR68HRbwL4/GiFQh/AjwL4esvH1DVfB/Dl0ddfBvArudd/bLQ64YsATnJpVqpoNE/i5wD8gZTy7+Z+xOvQECHE3VHGCEKIFQB/Bldzv34NwA+Pfm3yGqhr88MAflVyg7hapJR/XUq5J6V8B1f1/q9KKf8ieA0aJYRYE0LcUF8D+LMAfg+sj9rfBFII8YO4GnsOAfy8lPJvt3pAHhNC/DMAfxJXT1k+APA3Afw/AH4JwNsAHgL4ESnlp6NG/B/ganXbOYAfl1J+q4XD9ooQ4nsB/HsAv4vP5lr8DVzNO+J1aIAQ4r/D1STTEFcdxF+SUv7vQoj3cJXFeAvAbwH4n6WUAyHEMoB/jKv5YZ8C+FEp5YftHL1/RsNq/4uU8s/xGjRrdL7/1ejbCMA/lVL+bSHEbXS8Pmo9OCIiIiKySdvDakRERERWYXBERERElMPgiIiIiCiHwRERERFRDoMjIiIiohwGR0REREQ5DI6IiIiIchgcEREREeX8/3FI8MM295uIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(F,cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c4ffc26-a9f0-408c-80eb-723329c3b8ee",
   "metadata": {},
   "source": [
    "## 6. Solve and show"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "975aa11a-a4b7-49fb-869b-39843c1242e6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.9081958235744878e-17\n",
      "-1.9081958235744878e-17\n",
      "1.0983288221272578e-12\n"
     ]
    }
   ],
   "source": [
    "s_eff = np.linalg.inv(F) @ d\n",
    "s = A @ s_eff + s_N\n",
    "print(s_eff[-1])\n",
    "print(d[-1])\n",
    "print( np.linalg.norm(F @ s_eff - d))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8e4ab57-f7ab-4b05-a27f-61405d7d3079",
   "metadata": {},
   "source": [
    "## 7. Construct the flux field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "d4a4a3c5-081a-43ff-9f2a-1603428996bf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.19179212 -0.42309907]\n",
      "[ 0.19179212 -0.42309907]\n",
      "[ 0.15888522 -0.86499326]\n",
      "[ 0.15888522 -0.86499326]\n",
      "[ 0.25186634 -0.70030174]\n",
      "[ 0.25186634 -0.70030174]\n",
      "[ 0.17141968 -0.76678129]\n",
      "[ 0.17141968 -0.76678129]\n",
      "[ 0.24052252 -0.79747005]\n",
      "[ 0.24052252 -0.79747005]\n",
      "[ 0.19640065 -0.42217674]\n",
      "[ 0.19640065 -0.42217674]\n",
      "[ 0.22988955 -0.34988341]\n",
      "[ 0.22988955 -0.34988341]\n",
      "[ 0.27129175 -0.50776875]\n",
      "[ 0.27129175 -0.50776875]\n",
      "[ 0.1941863  -0.57996435]\n",
      "[ 0.1941863  -0.57996435]\n",
      "[ 0.15881071 -0.2       ]\n",
      "[ 0.15881071 -0.2       ]\n",
      "[ 0.26158548 -0.6045585 ]\n",
      "[ 0.26158548 -0.6045585 ]\n",
      "[ 0.2470484  -0.91249515]\n",
      "[ 0.2470484  -0.91249515]\n",
      "[ 0.20035552 -0.27204357]\n",
      "[ 0.20035552 -0.27204357]\n",
      "[ 0.19381801 -0.27196507]\n",
      "[ 0.19381801 -0.27196507]\n",
      "[ 0.18175402 -0.67034098]\n",
      "[ 0.18175402 -0.67034098]\n",
      "[ 0.22567741 -0.4987295 ]\n",
      "[ 0.22567741 -0.4987295 ]\n",
      "[ 0.18677721 -0.93015757]\n",
      "[ 0.18677721 -0.93015757]\n",
      "[ 0.24096258 -0.34880898]\n",
      "[ 0.24096258 -0.34880898]\n",
      "[ 0.25115496 -1.00682554]\n",
      "[ 0.25115496 -1.00682554]\n",
      "<visvis.wibjects.colorWibjects.Colorbar object at 0x0000019C957D4BB0>\n",
      "ERROR calling '_OnAxesPositionChange':\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\wibjects\\colorWibjects.py\", line 622, in _OnAxesPositionChange\n",
      "    self.position.x = axes.position.width+5\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\misc.py\", line 217, in fsetWithDraw\n",
      "    fset(self, *args)\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 1126, in fset\n",
      "    self._Update()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 928, in _Update\n",
      "    self._CalculateInPixels()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 997, in _CalculateInPixels\n",
      "    raise Exception(\"Can only calculate the position in pixels\"+\n",
      "Exception: Can only calculate the position in pixels if the owner has a parent!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "q_h = np.zeros((E,2)) \n",
    "\n",
    "for i, e in enumerate(element):\n",
    "    x1,y1 = node[e][0]\n",
    "    x2,y2 = node[e][1]\n",
    "    x3,y3 = node[e][2]\n",
    "    \n",
    "    A_e = 1/2*np.linalg.det(np.c_[node[e], np.ones((3,1))])\n",
    "    \n",
    "    Phi_e = -1/(2*A_e)*np.array([[x2-x3,x3-x1,x1-x2],\n",
    "                                [y2-y3,y3-y1,y1-y2]])\n",
    "    \n",
    "    if i in Eh:\n",
    "        L1 = np.sqrt((x2-x1)**2+(y2-y1)**2)\n",
    "        L2 = np.sqrt((x3-x2)**2+(y3-y2)**2)\n",
    "\n",
    "        left = np.linalg.inv(np.array([[(y1-y2)/L1,(x2-x1)/L1],\n",
    "                              [(y2-y3)/L2,(x3-x2)/L2]]))\n",
    "\n",
    "        Phi_L = np.c_[left, np.zeros((2,1))]\n",
    "        \n",
    "        print(Phi_e @ s[e])\n",
    "        Phi_e = np.c_[Phi_e, Phi_L @ Eh[i]]\n",
    "        \n",
    "        ey = np.append(e, -1)\n",
    "\n",
    "        q_e = Phi_e @ s[ey]\n",
    "        print(Phi_e @ s[ey])\n",
    "\n",
    "    else:\n",
    "        q_e = Phi_e @ s[e]\n",
    "    # q_e = Phi_e @ s[e]    \n",
    "    q_h[i] = q_e\n",
    "\n",
    "cfv.figure()\n",
    "\n",
    "# Draw the mesh.\n",
    "\n",
    "cfv.drawElementValues(\n",
    "    ev=q_h[:,1],\n",
    "    coords=coords,\n",
    "    edof=edof,\n",
    "    dofs_per_node=mesh.dofsPerNode,\n",
    "    el_type=mesh.elType,\n",
    "    title=\"Example 01\",\n",
    "    draw_elements=True\n",
    "        )\n",
    "cfv.showAndWait()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d7b83daa-f7c6-4f98-aeb6-0256c960d087",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.03504291382534708\n"
     ]
    }
   ],
   "source": [
    "fp1 = lambda x: -np.cosh(0.4)/np.sin(1)/np.sinh(1)*np.cos(x)\n",
    "fp2 = lambda y: -np.cos(0.6)/np.sin(1)/np.sinh(1)*np.cosh(y)\n",
    "fp3 = lambda x: np.cosh(0.6)/np.sin(1)/np.sinh(1)*np.cos(x)\n",
    "fp4 = lambda y: np.cos(0.4)/np.sin(1)/np.sinh(1)*np.cosh(y)\n",
    "\n",
    "\n",
    "print((fp1(0.6)-fp1(0.4))+(fp2(0.6)-fp2(0.4))+(fp3(0.4)-fp3(0.6))+(fp4(0.4)-fp4(0.6)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d717b115-2bdb-4587-a873-e504a0a5460d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1102230246251565e-16\n"
     ]
    }
   ],
   "source": [
    "fp1 = lambda x: np.cosh(0)/np.sin(1)/np.sinh(1)*np.cos(x)\n",
    "fp2 = lambda y: np.cos(1)/np.sin(1)/np.sinh(1)*np.cosh(y)\n",
    "fp3 = lambda x: -np.cosh(1)/np.sin(1)/np.sinh(1)*np.cos(x)\n",
    "fp4 = lambda y: -np.cos(0)/np.sin(1)/np.sinh(1)*np.cosh(y)\n",
    "\n",
    "\n",
    "print((fp1(1)-fp1(0))+(fp2(1)-fp2(0))+(fp3(1)-fp3(0))+(fp4(1)-fp4(0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "39f58a40-194f-4a63-98ba-75f206e02ab3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : GMSH -> Python-module\n",
      "Element=487, ae=[[ 24.99999994 -24.75316204   0.        ]], next <4,113>\n",
      "Element=664, ae=[[ 24.75316204   0.         -24.05371928]], next <140,113>\n",
      "Element=715, ae=[[ 24.05371928 -21.80623043   0.        ]], next <140,329>\n",
      "Element=929, ae=[[ 21.80623043   0.         -30.13294853]], next <584,329>\n",
      "Element=172, ae=[[-22.39270624  30.13294853   0.        ]], next <516,329>\n",
      "Element=153, ae=[[  0.          22.39270624 -26.36770403]], next <516,379>\n",
      "Element=21, ae=[[  0.          26.36770403 -22.50984033]], next <516,155>\n",
      "Element=810, ae=[[-21.13756962  22.50984033   0.        ]], next <149,155>\n",
      "Element=380, ae=[[  0.         -19.68341272  21.13756962]], next <154,155>\n",
      "Element=448, ae=[[ 19.68341272 -22.43218153   0.        ]], next <154,156>\n",
      "Element=582, ae=[[ 22.43218153   0.         -21.52188625]], next <158,156>\n",
      "Element=728, ae=[[ 21.52188625 -22.23017176   0.        ]], next <158,463>\n",
      "Element=102, ae=[[  0.         -21.10014342  22.23017176]], next <382,463>\n",
      "Element=109, ae=[[-21.930565     0.          21.10014342]], next <382,367>\n",
      "Element=142, ae=[[-20.50799238  21.930565     0.        ]], next <368,367>\n",
      "Element=40, ae=[[  0.          20.50799238 -19.71155648]], next <368,123>\n",
      "Element=612, ae=[[ 19.71155648   0.         -19.38177625]], next <552,123>\n",
      "Element=839, ae=[[  0.          19.38177625 -25.72221567]], next <552,15>\n",
      "Element=1016, ae=[[-19.99999995  25.72221567   0.        ]], next <14,15>\n",
      "21 [155 379 516 596]\n",
      "40 [123 367 368 596]\n",
      "102 [158 382 463 596]\n",
      "109 [382 367 463 596]\n",
      "142 [368 367 382 596]\n",
      "153 [379 329 516 596]\n",
      "172 [516 329 584 596]\n",
      "380 [149 154 155 596]\n",
      "448 [155 154 156 596]\n",
      "487 [ 84   4 113 596]\n",
      "582 [156 154 158 596]\n",
      "612 [123 368 552 596]\n",
      "664 [113   4 140 596]\n",
      "715 [113 140 329 596]\n",
      "728 [156 158 463 596]\n",
      "810 [149 155 516 596]\n",
      "839 [ 15 123 552 596]\n",
      "929 [329 140 584 596]\n",
      "1016 [ 14  15 552 596]\n",
      "size: 597 x 597\n",
      "rank: 596\n",
      "[[98], [1]] [[97, 96, 7, 95, 94, 93, 92, 6, 91, 90, 89, 88, 5, 87, 86, 85, 84, 4, 99], [27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 2]]\n",
      "first Neumann edge: [0.4  0.48] [0.4  0.52]\n",
      "constraint: [(98, 97), (97, 96), (96, 7), (7, 95), (95, 94), (94, 93), (93, 92), (92, 6), (6, 91), (91, 90), (90, 89), (89, 88), (88, 5), (5, 87), (87, 86), (86, 85), (85, 84), (84, 4), (4, 99)]\n",
      "constraint: (84, 4) [0.44 0.4 ] [0.4 0.4]\n",
      "first Neumann edge: [1. 0.] [1.   0.05]\n",
      "constraint: [(1, 27), (27, 28), (28, 29), (29, 30), (30, 31), (31, 32), (32, 33), (33, 34), (34, 35), (35, 36), (36, 37), (37, 38), (38, 39), (39, 40), (40, 41), (41, 42), (42, 43), (43, 44), (44, 45), (45, 2)]\n",
      "<visvis.wibjects.colorWibjects.Colorbar object at 0x0000017CE5B7AA30>\n",
      "ERROR calling '_OnAxesPositionChange':\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\wibjects\\colorWibjects.py\", line 622, in _OnAxesPositionChange\n",
      "    self.position.x = axes.position.width+5\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\misc.py\", line 217, in fsetWithDraw\n",
      "    fset(self, *args)\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 1126, in fset\n",
      "    self._Update()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 928, in _Update\n",
      "    self._CalculateInPixels()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 997, in _CalculateInPixels\n",
      "    raise Exception(\"Can only calculate the position in pixels\"+\n",
      "Exception: Can only calculate the position in pixels if the owner has a parent!\n",
      "\n",
      "True\n",
      "(array([  0,   0,   0, ..., 596, 596, 596], dtype=int64), array([  0,   8,  83, ..., 367, 368, 382], dtype=int64))\n"
     ]
    }
   ],
   "source": [
    "from fem import Poisson2D\n",
    "\n",
    "def create_geometry(el_on_curve=20):\n",
    "    # externel boundary ID\n",
    "    dir_id_1 = 100\n",
    "    dir_id_2 = 200\n",
    "    dir_id_3 = 300\n",
    "    neu_id_1 = 400\n",
    "    # hole boundaryd ID\n",
    "    hole_neu_id_1 = 500\n",
    "    hole_neu_id_2 = 600\n",
    "    hole_neu_id_3 = 700\n",
    "    hole_neu_id_4 = 800\n",
    "\n",
    "    # Geometry definition\n",
    "    g = cfg.Geometry()\n",
    "    g.point([0.0, 0.0], ID=0)\n",
    "    g.point([1.0, 0.0], ID=1)\n",
    "    g.point([1.0, 1.0], ID=2)\n",
    "    g.point([0.0, 1.0], ID=3)\n",
    "    g.point([0.4, 0.4], ID=4)\n",
    "    g.point([0.6, 0.4], ID=5)\n",
    "    g.point([0.6, 0.6], ID=6)\n",
    "    g.point([0.4, 0.6], ID=7)\n",
    "    g.spline([0, 1], el_on_curve=el_on_curve, marker=100, ID=0)\n",
    "    g.spline([1, 2], el_on_curve=el_on_curve, marker=101, ID=1) \n",
    "    g.spline([2, 3], el_on_curve=el_on_curve, marker=200, ID=2) \n",
    "    g.spline([3, 0], el_on_curve=el_on_curve, marker=300, ID=3)\n",
    "    g.spline([4, 5], el_on_curve=el_on_curve/4, marker=201, ID=4)\n",
    "    g.spline([5, 6], el_on_curve=el_on_curve/4, marker=301, ID=5)\n",
    "    g.spline([6, 7], el_on_curve=el_on_curve/4, marker=401, ID=6)\n",
    "    g.spline([7, 4], el_on_curve=el_on_curve/4, marker=501, ID=7)\n",
    "    g.surface([0, 1, 2, 3],[[4,5,6,7]])\n",
    "    cfv.drawGeometry(g)\n",
    "    cfv.showAndWait()\n",
    "\n",
    "\n",
    "    bcs = {100:[0],\n",
    "           200:[2],\n",
    "           300:[3],\n",
    "           101:[1],\n",
    "           201:[4],\n",
    "           301:[5],\n",
    "           401:[6],\n",
    "           501:[7],\n",
    "           }\n",
    "\n",
    "    for marker in bcs:\n",
    "        for i in bcs[marker]:\n",
    "            g.curve_marker(ID=i, marker=marker)\n",
    "    cfv.drawGeometry(g)\n",
    "    cfv.showAndWait()\n",
    "    return g\n",
    "\n",
    "\n",
    "def dir100(x,y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def dir200(x,y):\n",
    "    return np.sin(x)/np.sin(1)\n",
    "\n",
    "def dir300(x,y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def neu101(x,y):\n",
    "    return (1/np.tan(1))*np.sinh(y)/np.sinh(1)\n",
    "\n",
    "def neu201(x,y):\n",
    "    return (np.sin(x)/np.sin(1))*(np.cosh(0.4)/np.sinh(1))\n",
    "\n",
    "def neu301(x,y):\n",
    "    return -(np.cos(0.6)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "def neu401(x,y):\n",
    "    return -(np.sin(x)/np.sin(1))*(np.cosh(0.6)/np.sinh(1))\n",
    "\n",
    "def neu501(x,y):\n",
    "    return (np.cos(0.4)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "\n",
    "geometry = create_geometry(el_on_curve=20)\n",
    "Dirichlet = [100,200,300]\n",
    "Neumann = [101,201,301,401,501]\n",
    "hole = [[201,301,401,501]]\n",
    "\n",
    "BCfunc = {100:dir100, 200:dir200, 300:dir300,\n",
    "          101:neu101, 201:neu201, 301:neu301, 401:neu401,501:neu501}\n",
    "\n",
    "\n",
    "model = Poisson2D(geometry, Dirichlet, Neumann, hole, BCfunc)\n",
    "model.mesh(showMesh=True)\n",
    "model.efem()\n",
    "print(np.allclose(model.FU,FU))\n",
    "print(np.where(model.FU!=FU))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "79ff29a5-588a-4bea-9fa7-1212d5ae0652",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1598,)\n",
      "(1598,)\n"
     ]
    }
   ],
   "source": [
    "pos = np.where(model.FU!=FFF)\n",
    "print(pos[0].shape)\n",
    "print(pos[1].shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ee57bc83-d7e9-421b-a6a4-c959e11188bb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{(84,\n",
       "  4): array([[596],\n",
       "        [ -1]]),\n",
       " (14,\n",
       "  15): array([[596],\n",
       "        [  1]])}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.crossEdge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "78208e1d-f185-46e9-8116-3eba05f031b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from mayavi import mlab\n",
    "\n",
    "\n",
    "x = model.Node[:,0]\n",
    "y = model.Node[:,1]\n",
    "z = np.zeros_like(x)\n",
    "triangles1 = model.Element\n",
    "triangles2 = model.Element[list(model.crossElement.keys()),:]\n",
    "\n",
    "mlab.triangular_mesh(x, y, z, triangles1, color=(1,1,1),representation='wireframe' )\n",
    "mlab.triangular_mesh(x, y, z, triangles2, color=(1,0,0),representation='surface' )\n",
    "\n",
    "for cluster in model.NeuCluster:\n",
    "    i,_ = cluster[0]\n",
    "    x = model.Node[i][0]\n",
    "    y = model.Node[i][1]\n",
    "    z = 0    \n",
    "    mlab.points3d(x,y,z,1, color=(0,1,0),scale_factor=0.01)\n",
    "    for i, j in cluster:\n",
    "        x = model.Node[i][0]\n",
    "        y = model.Node[i][1]\n",
    "        z = 0\n",
    "        u = model.Node[j][0] - x\n",
    "        v = model.Node[j][1] - y\n",
    "        w = 0\n",
    "        mlab.quiver3d(x, y, z, u, v, w, scale_factor=1)\n",
    "\n",
    "for edge in model.crossEdge:\n",
    "    i,j = edge\n",
    "    x = [model.Node[i][0], model.Node[j][0]]\n",
    "    y = [model.Node[i][1], model.Node[j][1]]\n",
    "    z = [0,0]\n",
    "    mlab.plot3d(x, y, z, color=(0,0,1),tube_radius=0.005)\n",
    "        \n",
    "        \n",
    "mlab.view(0,0)\n",
    "# print(mlab.view())\n",
    "mlab.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "c920020e-b543-43b4-bf2b-07d819e1675e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "print(np.allclose(model.dU, dU))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "1e789dd1-81b1-401d-a221-891c892ff757",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "print(np.allclose(model.FU, FFF))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "554ae08f-ae06-4f49-872a-0f3b1841b6e0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys([(1, 27), (27, 28), (28, 29), (29, 30), (30, 31), (31, 32), (32, 33), (33, 34), (34, 35), (35, 36), (36, 37), (37, 38), (38, 39), (39, 40), (40, 41), (41, 42), (42, 43), (43, 44), (44, 45), (45, 2), (84, 4), (85, 84), (86, 85), (87, 86), (5, 87), (88, 5), (89, 88), (90, 89), (91, 90), (6, 91), (92, 6), (93, 92), (94, 93), (95, 94), (7, 95), (96, 7), (97, 96), (98, 97), (99, 98), (4, 99)])\n"
     ]
    }
   ],
   "source": [
    "print(model.Neumann.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "72e3da08-5193-40f5-afc6-6f92a6c237ba",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 97\n",
      "0 96\n",
      "0 7\n",
      "0 95\n",
      "0 94\n",
      "0 93\n",
      "0 92\n",
      "0 6\n",
      "0 91\n",
      "0 90\n",
      "0 89\n",
      "0 88\n",
      "0 5\n",
      "0 87\n",
      "0 86\n",
      "0 85\n",
      "0 84\n",
      "0 4\n",
      "0 99\n",
      "0 98\n",
      "1 27\n",
      "1 28\n",
      "1 29\n",
      "1 30\n",
      "1 31\n",
      "1 32\n",
      "1 33\n",
      "1 34\n",
      "1 35\n",
      "1 36\n",
      "1 37\n",
      "1 38\n",
      "1 39\n",
      "1 40\n",
      "1 41\n",
      "1 42\n",
      "1 43\n",
      "1 44\n",
      "1 45\n",
      "1 2\n"
     ]
    }
   ],
   "source": [
    "for i, cluster in enumerate(model.NeuCluster):\n",
    "    for e in cluster:\n",
    "        print(i,e[1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4233e0a1-b6e8-4467-9dc4-fce4a345a8ca",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EFEM",
   "language": "python",
   "name": "efem"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
